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Objectives. 

The  objectives  of  this  project  were  (1)  to  prove  convergence  theorems  for  probability-one 
homotopy  methods  applied  to  and  combined  optimal  model  order  reduction  and 

controller  synthesis  problems,  (2)  to  develop  a  robust,  fixed-structure  MATLAB  toolbox,  (3)  and 
to  extend  HOMPACK  to  deal  with  bifurcation  curve  tracking. 

Accomplishments/new  findings. 

For  three  different  formulations  of  the  optimal  model  order  reduction  problem  (optimal 
projection  equations,  input  normal  form  parametrization,  and  Ly  form  parametrization),  conver¬ 
gence  theorems  for  globally  convergent  probability-one  homotopy  algorithms  have  been  proved. 
Several  counterexamples  were  also  developed,  showing  that  the  results  are  sharp.  These  results 
complete  convergence  theory  for  optimal  model  order  reduction  homotopies.  Some  progress 
toward  homotopy  convergence  theory  for  combined  model  order  reduction  and  controller 

synthesis  was  also  made.  This  work  is  contained  in  an  Autom&tica  paper  under  review  and  in  Yuan 
Wang’s  Ph.D.  thesis.  The  current  version  of  the  Automatica  paper  is  attached  to  this  report. 

The  robust,  fixed-structure  MATLAB  toolbox  can  be  used  to  synthesize  fixed-structure  con¬ 
trollers  that  are  optimal  with  respect  to  a  given  performance  measure,  and  at  the  same  time  satisfy 
stability  and  robustness  constraints.  The  toolbox  can  handle  centralized  or  decentralized  compen¬ 
sators.  reduced  order  compensators,  or  compensators  with  repeated  gains,  all  in  a  common  format. 
The  available  performance  criteria  include  combined  maximum  entropy,  and  Popov. 

The  toolbox  has  been  tested  on  SUN,  DEC,  HP,  SGI,  and  IBM  UNIX  workstations,  and  UNIX 
viakc  files  are  provided  for  instaJlation  on  all  of  these  systems.  Documentation  for  the  toolbox  is 
appended  to  this  report.  Both  the  toolbox  and  documentation  are  available  at  the  URL: 

http://www.c3.vt.edu/  Itw/toolbox/ 

Personnel  supported. 

Computer  Science  M.S.  student  Kelly  O’Brien,  Ph.D.  student  Yuan  Wang,  and  the  PI  Layne 
Watson  were  supported  by  the  grant.  Faculty  and  students  associated  with  the  grant  include 
Dennis  Bernstein  (Michigan),  Scott  Erwin  (Michigan),  Yuzhen  Ge  (Butler),  and  Emmanuel  Collins 
(Florida  Ai:M). 
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tions  of  the  L*  optimal  model  order  reduction  problem”,  J.  Comput.  Appl.  Math.,  69  (1996) 
215-241. 

Y.  Ge,  L.  T.  Watson,  E.  G.  Collins,  Jr.,  and  D.  S.  Bernstein,  “Globally  convergent  homotopy  algo¬ 
rithms  for  the  combined  model  reduction  problem”,  J.  Math.  Systems,  Estimation, 

Control,  7  (1997)  129-155. 

Y.  Ge,  L.  T.  Watson,  E.  G.  Collins,  Jr.,  and  D.  S.  Bernstein,  “Probability-one  homotopy  algorithms 
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(1996)  187-208. 
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B.  B.  Lowekamp,  L.  T.  Watson,  and  M.  S.  Cramer,  “The  cellular  automata  paradigm  for  the 
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tracking”,  ACM  Trans.  Math.  Software,  22  (1996)  281-287. 

S.  Burgee,  A.  A.  Giunta,  V.  Balabanov,  B.  Grossman,  W.  H.  Mason,  R.  Narducci,  R.  T.  Haftka,  and 
L.  T.  Watson,  “A  coarse  grained  parallel  variable-complexity  multidisciplinary  optimization 
paradigm”,  Internat.  J.  Supercomputer  Appl.  High  Performance  Comput.,  10  (1996)  269—299. 

Gc.  L.  T.  Watson,  and  E.  G.  Collins,  Jr.,  “Cost-effective  parallel  processing  for  H^/H°°  con¬ 
troller  synthesis”,  Internat.  J.  Systems  Sci,  to  appear. 

M.  S.  Cramer,  S.  H.  Park,  and  L.  T.  Watson,  “Numerical  verificai,tion  of  scaling  laws  for  shock¬ 
boundary  layer  interactions  in  arbitrary  gases”,  J.  Fluids  Engrg.,  119  (1997)  67-73. 

A.  P.  Morgan,  L.  T.  Watson,  and  R.  A.  Young,  “A  Gaussian  derivative  based  version  of  JPEG  for 
image  compression  and  decompression”,  IEEE  Trans.  Image  Processing,  submitted. 

J.  F.  Monaco,  M.  S.  Cramer,  and  L.  T.  Watson,  “Supersonic  flows  of  dense  gases  in  cascade 
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E.  C.  Collins,  Jr.,  W.  M.  Haddad,  L.  T.  Watson,  and  D.  Sadhukhan,  “Probability-one  homotopy 
algorithms  for  robust  controller  synthesis  with  fixed-structure  multipliers”,  Internat.  J.  Robust 
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Math.  Software,  to  appear. 
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Convergence  Theory  of  Probability-one  Homotopies 
for  Model  Order  Reduction* 

Y.  WANGt,  D.  S.  BERNSTEIN*,  and  L.  T.  WATSON+ 

Theory  for  the  global  convergence  of  some  probability-one  homotopies  for  the  optimal 
model  order  reduction  problem  is  developed. 

Key  Words — Embedding;  globally  convergent  homotopy;  optimal  model  order  reduction;  probability-one  homotopy. 


Abstract — The  optimal  model  reduction  problem 

is  an  inherently  nonconvex  problem  and  thus 
provides  a  nontrivial  computational  challenge.  This 
paper  systematically  examines  the  requirements  of 
probability-one  homotopy  methods  to  guarantee  global 
convergence.  Homotopy  algorithms  for  nonlinear  systems 
of  equations  construct  a  continuous  family  of  systems 
and  solve  the  given  system  by  tracking  the  continuous 
curve  of  solutions  to  the  family.  The  main  emphasis 
IS  on  guaranteeing  transversality  for  several  homotopy 
mapft  based  upon  the  pseudogramian  formulation  of 
the  optimal  projection  equations  and  vernations  based 
upon  canonical  forms.  These  results  arc  essential  to 
the  probability-one  homotopy  approach  by  guarsinteeing 
good  numerical  properties  in  the  computational 
implementation  of  the  homotopy  algorithms. 

1.  INTRODUCTION 

Numerous  techniques  have  been  developed  to 
address  the  model  order  reduction  problem  with 
bolli  //•  and  criteria.  Model  reduction  from 
an  //*  perspective  is  considered  in  (Hyland  and 
Bornsiein,  1985)  and  (Beiratchart  et  al.,  1991). 
Baianrecl  truncation  and  associated  Hankel  norm 
reduction  theory  are  widely  used  in  practice 
to  provide  // '*'-suboptimad  solutions  ({Moore, 
1981),  (Glover,  1984),  (Zhou,  1995),  (Kabamba, 
1985b)),  More  recently,  convex  optimization 
inetliods  have  been  employed  iteratively  to 
approximate  solutions  to  the  nonconvex  problem 
(Grigoriadis,  1995).  These  techniques  are 
inherently  attractive  since  they  rely  only  upon 
convexity- based  procedures.  A  more  direct  albeit 
computationally  challenging  approach  is  to 
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apply  fixed-structure  optimization  (Hyland  and 
Bernstein,  1985),  (Haddad  and  Bernstein,  1989). 
Special  purpose  computational  methods  based 
upon  homotopy  techniques  have  been  developed 
for  this  problem  in  (Zigic  et  a].,  1993b),  (Ge  et 
a/.,  1994).  The  essential  difficulties  of  the  model 
reduction  problem  are  of  significant  interest  since 
techniques  developed  for  model  reduction  find 
immediate  application  to  the  closely  related 
problem  of  reduced-order  controller  synthesis 
((Hyland  and  Bernstein,  1984),  (Haddad  and 
Bernstein,  1990)). 

The  present  paper  is  concerned  with  the 
application  of  homotopy  methods  to  optimal 
model  reduction.  In  computational  practice, 
homotopy  methods  are  widely  used  for  nonconvex 
optimization  ((Watson,  1990),  (Watson  and 
Haftka,  1989)).  Homotopy  methods,  in 
particular,  probability-one  homotopy  methods, 
have  global  convergence  properties  that  are  often 
advantageous  in  comparison  to  locally  convergent 
methods  such  as  quasi-Newton  methods  ((Chow 
et  ai.,  1978),  (Watson,  1989),  (Watson,  1986)). 
Under  suitable  hypotheses,  probability-one 
homotopy  methods  are  guaranteed  to  converge 
globally  (from  an  arbitrary  starting  point)  to 
a  solution  of  a  nonlinear  system  of  equations. 
The  nomenclature  ‘‘probability-one”  is  well 
established  in  the  mathematical  literature  and 
reflects  the  generic,  measure  theoretic  properties 
of  the  algorithms  rather  than  stochastic  aspects. 

The  goal  of  the  present  paper  is  to 
systematically  examine  the  requirements  of 
probability-one  homotopy  methods  to  guarantee 
global  convergence.  The  crucial  requirements 
are  1)  transversality  and  2)  boundedness.  As 
discussed  in  Section  2,  transversality  implies  the 
existence  of  eind  the  ability  to  track  a  zero 
curve  of  the  homotopy  map,  while  boundedness 
is  equivalent  to  the  existence  of  solutions  to 
the  model  reduction  problem.  The  existence  of 
optimal  reduced-order  models  follows  from 
the  results  in  (Spanos  et  af.,  1990).  The  main 
emphasis  in  the  present  paper  is  on  guaranteeing 
transversality  for  several  homotopy  maps  based 
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upon  the  pseudogreunian  formulation  of  the 
optimal  projection  equations  and  specialized 
formulations  based  upon  canonical  forms.  These 
results  are  essential  to  the  probability-one 
homotopy  approach  by  guaranteeing  good 
numerical  properties  (explained  in  (Watson  et 
ai.,  1987))  in  the  computational  implementation 
of  the  homotopy  algorithms.  Numerical 
comparisons  with  other  approaches  have  been 
done  elsewhere  ((Ge  et  ai.,  1996),  (Zigic  et  a/., 
1993b)),  and  are  not  the  objective  of  the  present 
paper. 

The  contents  of  the  paper  are  as  follows. 
After  stating  the  model  reduction  problem 
in  Section  2,  we  then  provide  a  brief 
review  of  probability-one  homotopy  theory  in 
Section  3.  The  transversality  assumption  of 
probability-one  homotopy  theory  is  then  proven 
in  Section  4  for  several  canonical  forms. 
Next,  it  is  shown  by  example  in  Section  5 
that  the  boundedness  assumption  required  by 
probability-one  homotopy  theory  is  not  aJways 
satisfied  by  the  pseudogramian  formulation  of 
the  optimal  projection  equations  and  by  some 
formulations  based  upon  canonical  forms.  Then 
it  IS  shown  that  for  a  reformulation  of  the 
pseudogramian  optimal  projection  equations  in 
complex  projective  space  using  homogeneous 
transformations,  the  boundedness  assumption 
holds  and  thus  convergence  of  the  homotopy 
algorithm  to  a  solution  (in  complex  projective 
space)  is  guaranteed.  The  numerical  results  in 
(Ge  et  a/.,  1996)  and  (Zigic  et  a/.,  1993b)  show 
that,  in  praettre.  it  is  not  necessary  to  track 
the  homotopy  zero  curves  in  complex  projective 
space.  Section  6  concludes. 

2.  //^  OPTIMAL  MODEL  ORDER  REDUCTION 

I’he  //^  optimal  model  order  reduction 
problem  can  be  formulated  as  follows:  given  the 
nth-order  asymptotically  stable,  controllable  and 
observable  linear  time- invariant  continuous- time 
system 

i(t)  =  Ax{t)  +  (2.1) 

y(l)  =  Cx{t),  (2.2) 

where  >4  6  R”*",  fl  6  R”*'"*,  and  C  G  R***", 
and  given  <  n,  find  an  Umth-order 

reduced-order  model 

im(i)  =  A^Xrr^it)  +  (2.3) 

l/m(0  ~  Cn,®m(0i  (2-^) 

where  Am  €  is  asymptotically  stable, 

Bm  €  Cm  €  R*’'"”,  which  minimizes 

the  quadratic  model-reduction  criterion 
( -^m  1  Bm  1  Cm  )  = 

lim  ^[(y(0  -  y„(<))^ii(y(<)  -  y„.(<))K2.5) 


where  the  input  u{t)  is  white  noise  with  positive 
definite  intensity  7,  and  ii  is  a  positive  definite 
weighting  matrix.  Throughout,  all  positive 
semidefinite  and  positive  definite  matrices  are 
assumed  to  be  symmetric. 

To  guarantee  that  J  is  finite,  a 
solution  {Am,  Bm,  Cm)  is  sought  in  the 
set  5  =  {(A  m  I  Bm ,  Cm)  •  Am  is  asymptotically 
stable,  {Am,  Bm)  is  controllable  and  {Am,  Cm)  is 
observable}.  In  this  case  the  quadratic  model 
reduction  criterion  (2.5)  is  given  by 

J{Am,  Bm,Cm)  ^  (^’6) 

where 


C  =  {C  ^Cm),  R  =  C^RC 


and 

g=:  /  e^^BVB'^e^'^Ui, 

Jo 

which  is  the  unique  solution  of  the  Lyapunov 
equation 

AQ  +  QA^  -h  BVB'^  =  0.  (2.7a) 

For  future  reference  define  P  by 

=  (2.76) 

eoid  partition  P,  Q  as 


in  conformance  with  A. 

The  following  theorems  and  lemmas  from 

(Haddad  and  Bernstein,  1989),  (Hyland  and 

Bernstein,  1985)  will  be  needed  in  Section  4. 

Lemma  2,1.  Let  positive  semidefinite 

QyP  £  satisfy 

rank{(3)  =  rank(P)  =  rank((5-P)  =  rim,  (2.8) 
where  rim  <  n.  Then  there  exist  nonsingular 
W  G  cind  positive  definite  diagonal 

E  G  such  that 

2) 

Lemma  2.2.  Let  positive  semidefinite 

Q,  PgR”^”  satisfy  the  rank  conditions  (2.8), 
where  rxm  <  n.  Then,  there  exist  G,  F  G 
and  positive  semisimple  Af  G  unique 

up  to  a  change  of  basis  in  R”*",  such  that 

QP  =  G'^MT,  rG’’  =  /„„.  (2.9) 

Theorem  2.3.  Suppose  {Am,  Bm,  Cm)  G  S  solves 
the  optimal  model-reduction  problem.  Then 
there  exist  positive  semidefinite  matrices  Q, 
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P  G  R"^"  satisfying  (2.8)  and  such  that  Am,  Bm 
and  Cm  are  given  by 

>1^  =  TAG'^,  Bm  =  rP,  Cm  =  CC^,  (2.10) 

and  such  that,  with  r  =  C^T,  the  following 
conditions  are  satisfied: 

t[AQaQA^  +  BVB'^]  =  ^,  (2.11) 

[A^  P  +  PAaC'^  RC\t  =  Q.  (2.12) 

3.  PROBABILITY-ONE  GLOBALLY  CONVERGENT 
HOMOTOPIES 

A  homotopy  is  a  continuous  map  from  the 
interval  [0,1]  into  a  function  sp2u:e,  where  the 
continuity  is  with  respect  to  the  topology  of 
the  function  space.  Intuitively,  a  homotopy  p(A) 
continuously  deforms  the  function  p(0)  =  g  into 
the  function  p(l)  =  /  as  A  goes  from  0  to  1.  In 
this  case,  /  and  g  are  said  to  be  homotopic. 
Homotopy  maps  are  fundamental  tools  in 
topology,  and  provide  a  powerful  mechanism  for 
defining  equivalence  classes  of  functions. 

Homotopies  provide  a  mathematical  formalism 
for  describing  an  old  procedure  in  numerical 
analysis,  variously  known  as  continuation, 
incremental  loading,  amd  embedding.  The 
continuation  procedure  for  solving  a  nonlinear 
system  of  equations  f{x)  =  0  stairts  with  a 
(generally  simpler)  problem  =  0  whose 

solution  xo  is  known.  The  continuation 
procedure  is  to  track  the  set  of  zeros  of 

p{X,x)  =  Xnx)+{l-X)g(x)  (3.1) 
as  A  is  increased  monotonically  from  0  to 
1,  starting  at  the  known  initial  point  (0,  xq) 
satisfying  p(0,Xo)  =  0.  Each  step  of  this  tracking 
process  is  done  by  starting  at  a  point  (A,  x)  on 
the  zero  set  of  p.  fixing  some  AA  >  0,  and  then 
solving  p(A-f  AA,x)  =  0  for  x  using  a  locally 
convergent  iterative  procedure,  which  requires  an 
invertible  Jacobian  matrix  £)xp(A -f  AA,x).  The 
process  stops  at  A  =  1,  since  /(x)  =  p{i,i)  =  0 
gives  a  zero  x  of  /(x).  Note  that  continuation 
assumes  that  the  zeros  of  p  connect  the  zero  xq 
of  ^  to  a  zero  x  of  /,  and  that  the  Jacobian 
matrix  Drp{A,  x)  is  invertible  along  the  zero  set 
of  p;  these  are  strong  assumptions,  which  are 
frequently  not  satisfied  in  practice. 

Continuation  can  fail  because  the  curve  7  of 
zeros  of  p(A,  x)  emanating  from  {0,xo)  may  (1) 
have  turning  points,  (2)  bifurcate,  (3)  fail  to 
exist  at  some  A  values,  or  (4)  wander  off  to 
infinity  without  reaching  A  =  1.  Turning  points 
and  bifurcation  correspond  to  singular  Dxp{X,x). 
Generalizations  of  continuation  known  as 
homotopy  methods  attempt  to  deal  with  cases  (1) 
and  (2),  and  allow  tracking  of  7  to  continue 
through  singularities.  In  particular,  continuation 


monotonically  increases  A,  whereas  homotopy 
methods  permit  A  to  both  increase  and  decrease 
along  7.  Homotopy  methods  can  also  fail  via 
cases  (3)  or  (4). 

The  map  p(A,x)  connects  the  functions  g{x) 
and  /(x),  hence  the  use  of  the  word  “homotopy.” 
In  general  the  homotopy  map  p(A,x)  need  not 
be  a  simple  convex  combination  of  g  and  / 
as  in  (3.1),  and  can  involve  A  nonlinearly. 
Sometimes  A  is  a  physical  parameter  in  the 
original  problem  /(x;  A)  =  0,  where  A  =  1  is  the 
(nondimensionalized)  value  of  interest,  although 
“artificial  parameter”  homotopies  are  generally 
more  computationally  efficient  than  “natural 
parameter”  homotopies  p(A,x)  =  /(x;A).  An 
example  of  an  artificieJ  parameter  homotopy  map 
is 

p{X,x)  =  Xf{x;X)  +  {l-X){x-a),  (3.2) 

which  satisfies  p(0,a)  =  0.  The  name  “artificial” 
reflects  the  fact  that  solutions  to  p(A,  x)  =  0  have 
no  physical  interpretation  for  A  <  1.  Note  that 
p(A,  x)  in  (3.2)  has  a  unique  zero  x  =  a  at  A  =  0, 
regardless  of  the  structure  of  /(x;  A). 

All  four  shortcomings  of  continuation  and 
homotopy  methods  have  been  overcome  by 
probability-one  homotopies,  proposed  in  1976 
by  Chow,  Mallet-Paret,  and  Yorke  (Chow  et 
ai.,  1978).  The  supporting  theory,  based  on 
differential  geometry,  will  be  reformulated  in  less 
technical  jargon  here. 

Definition  3.1.  Let  U  C  and  V  C  R^  be 
open  sets,  and  let  p  :  x  [0, 1)  x  F  —)*  R^  be  a 
map.  p  is  said  to  be  transversal  to  zero  if  the 
p  X  (m  4*  1  +  p)  Jacobian  matrix  Dp  has  full  rank 
on  p"^(0). 

The  requirement  is  technical,  and  part  of 
the  definition  of  transversality.  The  basis  for  the 
probability-one  homotopy  theory  is: 

Theorem  3.2  (Parametrized  Sard’s  Theorem) 
(Chow  et  a/.,  1978).  Let  p  :  x  [0, 1)  x  1/  ->  R^ 
be  a  map.  If  p  is  transversal  to  zero,  then  for 
almost  allaeU  the  map 

Pa{X,x)  =  p{a,X,x) 
is  also  transversal  to  zero. 

To  discuss  the  import  of  this  theorem,  take 
U  =  R”*,^  =  R^,  and  suppose  that  the 
map  p  :  R”'  x  [0, 1)  x  R^  -4  R^  is  transversal 
to  zero.  A  straightforward  application  of  the 
implicit  function  theorem  yields  that  for  almost 
all  a  E  R”',  the  zero  set  of  pa  consists  of 
smooth,  nonintersecting  curves  which  either  (1) 
are  closed  loops  lying  entirely  in  (0,l)xRP, 
(2)  have  both  endpoints  in  {0}  x  R^,  (3)  have 
both  endpoints  in  {1}  x  RP,  (4)  are  unbounded 
with  one  endpoint  in  either  {0}  x  R^  or  in 
{1}  X  RP,  or  (5)  have  one  endpoint  in  {0}  x  R^ 
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and  the  other  in  {1}  x  R^.  Furthermore,  for 
almost  all  a  E  R*”,  the  Jacobian  matrix  Dpa  has 
full  rank  at  every  point  in  The  goal  is 

to  construct  a  map  pa  whose  zero  set  has  an 
endpoint  in  {0}  x  RP,  and  which  rules  out  (2) 
and  (4).  Then  (5)  obtains,  and  a  zero  curve 
starting  at  (0,  xq)  is  guaranteed  to  reach  a  point 
(l,x).  All  of  this  holds  for  almost  all  a  G  R”*, 
and  hence  with  probability  one  (Chow  et  ai., 
1978).  Furthermore,  since  a  G  R”*  can  be  almost 
any  point  (and,  indirectly,  so  cam  the  starting 
point  xo),  an  algorithm  based  on  tracking  the 
zero  curve  in  (5)  is  legitimately  called  globally 
convergent.  This  discussion  is  summarized  in  the 
following  theorem. 

Theorem  3.3.  Let  f  '."RP  —¥  RP  be  a  map, 
p  :  R""  X  [0, 1)  X  RP  R^  a  map,  and 
Pa{X,x)  =  />(a,A,x).  Suppose  that 

( 1 )  p  is  transversal  to  zero, 
and,  for  eeuch  fixed  a  G  R'”, 

(2)  pa{0,x)  =  0  has  a  unique  solution  Xo, 

(:i)  Pa{i,x)  =  f(x)  (leRP). 

Then,  for  almost  all  a  G  R"',  there  exists  a  zero 
curve  7  of  Pa  emanating  from  (0,xo),  along 
which  the  Jacobian  matrix  Dpa  has  full  rank.  If, 
in  addition, 

(4)  p“'(0)  is  bounded, 

then  7  reaches  a  point  (l,x)  such  that  /(x)  =  0. 
Furthermore,  if  Df{x)  is  invertible,  then  7  has 
finite  arc  length. 

Any  algorithm  for  tracking  7  from  (0,Xo)  to 
(l,x),  based  on  a  homotopy  map  satisfying 
the  hypotheses  of  Theorem  3.3,  is  called 
a  globally  convergent  probability- one  homotopy 
algorithm.  Of  course  the  practical  numerical 
details  of  tracking  7  are  nontrivial,  and 
have  been  the  subject  of  twenty  years  of 
research  in  numerical  analysis.  Production 
quality  software  called  HOMPACK  (Watson 
et  ai.,  1987)  exists  for  tracking  7.  The 
distinctions  between  continuation,  homotopy 
methods,  and  probability-one  homotopy  methods 
are  subtle  but  worth  noting.  Only  the 
latter  are  provably  globaJly  convergent  and 
(by  construction)  expressly  avoid  dealing  with 
singularities  numerically,  unlike  continuation  and 
homotopy  methods  which  must  explicitly  handle 
singularities  numerically. 

The  purpose  of  this  paper  is  to  prove  or 
disprove  properties  (l)-(4)  of  Theorem  3.3  for 
some  homotopy  maps  that  have  been  proposed 
for  the  optimal  model  order  reduction 

problem,  and  which  have  been  successful  in 
practice.  Assumptions  (2)  and  (3)  in  Theorem 
3.3  are  usually  achieved  by  the  construction 
of  p  (such  as  (3.2)),  and  are  straightforward 
to  verify.  Although  assumption  (1)  is  trivial 
to  verify  for  some  maps,  for  the  model 


order  reduction  homotopies  the  verification  is 
nontrivial.  Assumption  (4)  is  typically  very  hard 
to  verify,  and  often  is  a  deep  result,  since  (l)-(4) 
holding  implies  the  existence  of  a  solution  to 
/(x)=0. 

Note  that  (1)”(4)  are  sufficient,  but  not 
necessary,  for  the  existence  of  a  solution  to 
f[x)  =  0,  which  is  why  homotopy  maps  not 
satisfying  the  hypotheses  of  Theorem  3.3  can 
still  be  very  successful  on  practical  problems.  If 
(l)-(3)  hold  and  a  solution  does  not  exist,  then 
(4)  must  fail,  and  nonexistence  is  manifested 
by  7  going  off  to  infinity.  Properties  (l)-(3) 
are  important  because  they  guarantee  good 
numerical  properties  along  the  zero  curve  7, 
which,  if  bounded,  results  in  a  globally  convergent 
algorithm.  If  7  is  unbounded,  then  either  the 
homotopy  approach  (with  this  particular  p)  has 
failed  or  /(x)  =  0  has  no  solution. 

4.  TRANSVERSALITY  OF  HOMOTOPIES  FOR 
OPTIMAL  MODEL  ORDER  REDUCTION 

This  section  proves  that  three  homotopies  /)(a, 
A,x)  which  have  been  used  in  (Zigic  et  ai., 
1993b)  and  (Ge  et  a/.,  1994)  for  the  optimal 
model  order  reduction  problem  are  transversal 
to  zero,  the  first  requirement  of  Theorem  3.3. 
An  overview  and  compa^’kon  of  these  homotopy 
maps  is  in  (Ge  et  a/.,  1996).  The  analysis 
concerns  (2.11)  and  (2.12)  where  Q  and  P  are 
positive  semidefinite  matrices  satisfying  (2.8). 

4.1.  Transversality  of  homotopies  based  on 
decompositions  of  pseudogram  ians 

Since  Q  and  P  satisfy  (2.8),  there  exists 
invertible  W  E  R”^"  and  positive  definite 
diagonal  E  G  such  that  (Hyland  and 

Bernstein,  1985) 

g  =  W'(o  J)  = 

where 

w=  {Wi  W2),  =  =  ([^^). 

Premultipling  (2.11)  by  Ui  and  postmultipiying 
(2.12)  by  Wi  yields  (recall  that 
T=:G^r=:WiUi) 

UiAWii:Wj  +  LW^A^  +  UiBVB'^=0,[A.l) 
A'^U‘[i:  +  Ufi:UiAWi  +  C'^  RCWi  =  0.{4.2) 
A  constraint  from  W~^  =  U  is 
UiWi  -1  =  0. 


(4.3) 
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The  matrix  equations  (4.1)-'(4.3)  contain 
2 nn,7i  +  scalar  equations.  However,  the 
only  unknowns  in  (4.1)-(4,3),  namely  Wi,  i7i, 
and  diagonal  E,  contain  2nnm-\-Tim  variables. 
Hence,  some  other  formulation  is  necessary 
in  order  to  make  an  exact  match  between 
the  number  of  equations  and  the  number  of 
unknowns.  Following  (Zigic  et  aJ.,  1993b),  all 
elements  of  E  are  considered  as  unknowns,  giving 
the  same  number  of  equations  as  unknowns.  The 
structure  of  the  problem  is  such  that  E  will  turn 
out  to  be  symmetric,  so  it  can  be  diagonalized  to 
produce  the  decomposition  of  Q  and  P  described 
above. 

The  approax:h  in  (Zigic  et  aJ.,  1993b),  analyzed 
next,  uses  the  homotopy  map 

/>a(A,x)  =A/(i) +  (1-A)fli(x;a),  (4.4) 

where  the  initial  problem  pa(0,x)  =  <7(x;o)  =  0 
ha.s  an  easily  obtained  unique  solution  and  the 
final  problem  (4.1)-(4.3)  is  pa{lyx)  =  f{x)  =  0. 
/  and  g  are  displayed  in  (4.4)  simply  to  point 
out  that  the  map  pa(A,x)  can  be  viewed  as 
a  rt>Mvt*x  combination  of  two  other  ma|>s.  For 
notational  convenience  later  when  displaying 
Jarr>bian  matrices  the  order  of  the  variables  is 
limr*  forth  taken  as  A,  x,  a.  Lret 

A{X)^A,  B(X)  =  XBVB'^ +  {l-X)Bi, 

C(A)  =  AC^/ZC+(1-A)C<, 

where  ii,  =  B{0)  and  C,  =  C(0)  are  matrices 
defining  the  initial  problem  at  A  =  0,  and 

correspond  to  the  parameter  vector  a  in  Theorem 
3  3  Define 

(Fi(X,x,a)\ 

fuiX.T)  =  piX.x.a)  =  I  Fi(A,x,a) 
\F3(A,x,a)/ 

in  (‘l.dj  by 

FliX.x.a)  =  +  >1^(A) 

+  fi/i|A).  (4,5) 

F,(A,r,«)  =  A^{X)V:[  E  +  Uf  EF,  yl(A)  W'l 
+  C’(A)1V,.  (4  6) 

F3{X,x,a)  =  Ui\Vi-I,  (4.7) 

where 


is  the  generic  parameter  vector  in  Theorem  3.3 
and  in  (4.4), 

/y/ec(Wi)\ 
x=  Vec(t/,) 

V  Vec  (E)  ) 

denotes  the  independent  variables  Wi  €  R"*'"’", 
Ui  €  R"”'*'",  E  €  R""**'”"*  corresponding  to  x  in 
Theorem  3.3,  and  A,  B,C,V,  R  are  constants  as 
in  Section  2. 


The  Jacobian  matrix  of  p(A,x,a)  has 

2  n  rim  +  fim  2  n*  +  2  n  Tim  +  Om  +  1 

columns.  Rows  1  through  n  rim  correspond 
to  (4.5),  rows  nrim  +  l  through  2nnm 
correspond  to  (4.6),  and  rows  2nnm  +  1  through 
2nnm  +  nm  correspond  to  (4.7).  The  first 
column  corresponds  to  the  derivatives  with 
respect  to  A,  columns  2  through  n  + 1 
correspond  to  the  derivatives  with  respect  to  Wi , 
columns  nnm  +  2  through  2  n  +  1  correspond 
to  the  derivatives  with  respect  to  Ui,  columns 
2  nnm  +  2  through  2  n  rim  +  ^m  +  1  correspond 
to  the  derivatives  with  respect  to  E,  columns 
2  n  rim  +  "m  +  2  through  2  n  Um  +  +  1 

correspond  to  the  derivatives  with  respect 
to  Bi,  and  columns  2nnm  +  Hm  +  ”^  +  2 
through  2  n  rim  +  rim  +  2  +  1  correspond  to  the 

derivatives  with  respect  to  Cii 

Dp{X,x,a)  =  {Dxp  DwiP  DuiP 

Dzp  DsiP  Dc,p)-  (4.8) 
Since  F3(X,x,a)  does  not  depend  upon  A,  Bi, 
and  Ci,  it  follows  that 

DxF3{X,x,a)  =  0, 

■DBif^3(A,x,a)  =  0, 

Dc.Fsi^tXtO)  =  0, 


and  simil^a'ly 

DciFi{X,x,a)  =  DB,F2{X,x,a)  =  0. 

Thus 

Dp{X, X, a)  =  DpiX,  WuUu^, Bi,Ci) 

/DxFi  D^Fi  D,Fi' 


=  DxF2 
\  0 
(DxFi 

=  DxF2 
0 


DrF2 
Dr  Fa 


\ 


Dx,Fi 

Dr:F2 

D^Fa 


(4.9) 


The  following  lemma  will  be  used  in  the  proof 
of  Theorem  4.2. 


Lemma  4.1.  Let  X  €  R''"'’  and  A  €  R"'"'", 
B  e  R”***  be  differentiable  with  respect  to  x,j 
for  1  <  i  <  p,  1  <  j  <  q.  Then 


d 

dXij 


{AB)  =  ( 


d 

dxij 


and  for  constant  M,  interpreting  the  derivative 
Dx  as  Dvec(X)> 


Dx [MX)  =  I®  M,  Dx [XM)  =  ®  L 


The  proof  of  Lemma  4.1  is  straightforward 
calculus. 

Theorem  4.2.  The  homotopy  map  given  by 
(4.5)-(4.7)  is  transversal  to  zero  (for  0  <  A  <  1). 
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Proof.  To  prove  that  Dp{\x,a)  given  in  (4.9) 
has  full  rank,  i.e., 

rank  (I^p(A,x,a))  =  2nnm  + 
it  suffices  to  prove  that 

Ta.nk{DxFs)  =  xsink{DwiF3  Du^Fz  D^iFs) 

=  n^,  (4.10) 

rank(DaPi)  =  rank(f>B<-Pa  0)  =  (4-11) 

rank(f>ai^2)  =  rank(0  DciF2)  =  (4.12) 

The  meaning  of  expressions  like  D^Fz  is 

ambiguous  until  some  ordering  is  specified  for  the 
components  of  the  matrices  E  and  F3.  Hereafter, 
whichever  ordering  is  notationally  convenient  is 
used.  If  unspecified,  the  standard  ordering  by 
columns  (Vec)  is  assumed. 

Using  Lemma  4.1,  ordering  Ui  and  Fz  by  rows, 

Dr,F3{X,x,a)  =  Du.iUiWi)  =  /„„  ® 

(4.13) 

and  ordering  Wi  and  F3  by  columns, 

Ow,F3{X,x,a)  =  Dwi{^iW\)  —  In„  ®  Ui- 

(4.14) 

Siiirf  r'lH’i  =  I,  by  Sylvester’s  inequdity, 
rank  (Ui)  =  rank  (VTi)  —  Hni, 
and  tliercfore 

rank(DrF3)  =  rank(Z)(;jF3) 

=  rank(Dw,  F3)  =  n^, 
which  is  (4.10). 

I 'sing  Lemma  4.1,  ordering  Bi  and  Fj  by 
columns  yields 

Di,,Fi(X.x.a)  =  DbAUiB(X)) 

=  {l-  X)DbAUiB,) 

=  {\-X)I„®Uu  (4.15) 

and  using  (4.15)  for  A  <  1  yields 

rank{  Db,  Fi  )  =  nrim. 

Similarly,  ordering  Ci  and  Fj  by  rows, 
Dc.F2(X,x,a)  =  DcAC{X)Wi) 

=  {1  -  X)DcACiWi) 

=  (1  -  A)/„  ®  W^,  (4.16) 

so  for  A  <  1 

rank(Dc.F2)  =  nrim- 

This  completes  the  proof  of  (4.10)-(4.12),  and 
the  proof  that  the  homotopy  map  (4.5)“(4.7)  is 
transversal  to  zero  for  all  0  <  A  <  1.  Q.  E.  D. 


4.2.  Transversality  of  bomotopies  based 
on  input  normal  form 

The  following  theorem  from  (Kabamba,  1985a) 
is  needed  to  present  the  homotopy  method  for 
the  input  normal  form. 

Theorem  4,3,  Suppose  is 

asymptotically  stable  and  minimal.  Then  there 
exist  a  simileirity  transformation  U  and  a  positive 
definite  matrix  Q  =  diag(a;i,  •  •  ‘  ,a;n,„)  such  that 
Am  =  U^^AmU,  Bm  =  and  Cm  =  CmU 

satisfy 

Am+Am  +  BmVB'^  =  0, 

AlQ^QAm-\-Cf^RCm  =  0^ 

In  addition,  if  the  w,*  are  distinct, 

{Amh  =  -\{B,„VBl),,, 

{ClRCmh 

{Brr.VBr)^’ 

Wj  —  Wi 

(4.18) 


(>lm) 


Definition  4.3.1,  The  triple  {Am^Bm^Cm) 
satisfying  (4.17)  and  (4.18)  is  said  to  be  in  input 
normal  form. 

The  utility  of  the  input  normal  form 

(4.17) -(4.18)  lies  in  using  Bm  and  Cm  as  the 
independent  variables,  and  then  being  able  to 
recover  Am  uniquely  from  Bm  and  Cm-  The 
number  of  variables  in  Bm  and  Cm  is  nyn(m  +  /), 
the  minimum  number  of  variables  possible  to 
describe  any  reduced  order  model,  and  thus  the 
input  normal  form  parametrization  is  referred  to 
as  a  “minimal  parametrization.”  If  Wi  =  wj  for 
some  I  ^  jy  then,  regardless  of  (4.17)  holding, 

(4.18)  fails  to  permit  the  unique  recovery  of  Am- 
Under  the  assumption  that  the  solution  (^4^, 

Bmi  Cm)  being  sought  exists  in  input  normal 
form,  the  only  independent  variables  are  Bm  and 
Cm  I  and  in  this  case  the  domain  is 
{(Am^BmyCm)  ^  Am  is  asymptotically  Stable, 
{AmyBm^Cm)  is  minimal  and  in  input 
normal  form}. 

Now  for  {Am,  Bm,  Cm)  in  input  normal  form, 
the  cost  function  can  be  written  as 

J(A„.,B„„C„,)  =  tr  {QjRi),  (4.19) 


Remark  4.2.1.  One  can  use  more  variables  in  the 
parameter  vector  a,  e.g.,  i4(A)  =  AA  -h  (1  “  A)>1|, 
without  affecting  the  full  rank  properties. 


where  Qj  is  a  symmetric  and  positive  definite 
matrix  satisfying 

AiQi  +  QiAJ  ~\-Vj  =  0 


(4.20) 
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and 


Rj  = 


Vi  = 


(  CFRC 

(  BVB'^ 
\BmVB^ 


-Cf^RCr„\ 
ClRCm  )  ’ 
BVBl  \ 
BmVBlJ- 


(4.21) 


Qj  can  be  written  as 

where  Qi  €  R"’'",  Q12  €  and 

Q2  G 

Minimizing  (4.19)  under  the  constraints  (4.17) 
and  (4.20)  leads  to  the  Lagrangian 

L(A,n,  B,n,Cm,^,Ql,  Me,  Mo,  Pi) 

=  tr  +  {Am  +  Am  +  BmV  Bm)  Me 

+  {AlQ  +  QAm  +  ClRCm)Me 

+  {AiQi  +  QiAJ  +  K/)P/j, 


whf're  the  symmetric  matrices  Me,  and  Pi 
are  Lagrange  multipliers. 

Setting  dL/dQi  =  0  gives  an  equation  for  Pj 
sirmiar  to  (4.20)  for  P, 

Aj  P,  +  PiAi  +  Ri  =  0,  (4.23) 


where  P/  is  symmetric  positive  definite  and  can 
be  partitioned  similarly  to  Qj  as 

^'=ik  ft)' 

The  matrices  Afc  and  Afo  satisfy  (Davis  et  ai., 
1902) 

A/,  = -(is  +  fiM,),  (4.25) 

V"  "i/ii  j»i 


where  S  =  2(Pj^0i3  +  P2Q2)- 
Setting  dL/dBm  =  0  and  dL/dCm  =  0  gives 
2(P,^2B  +  P2Bm)V  +  2MeBmV  =  0.  (4.28) 

2R{CmQ2  -  CQ12)  +  2RCmMe  =  0.  (4.29) 

Observe  that  P/  through  (4.23)  and  Qi  through 
(4.20)  depend  on  Bm  and  Cm  as  does  Am 
through  (4.18).  Similarly  Me  through  (4.25) 
and  Mo  through  (4.26)-(4.27)  depend  on  Bm 
and  Cm-  Thus  everything  in  (4.28)“(4.29)  is  a 
function  of  Bm  and  Cm-  Use  the  homotopy  map 
structure  of  (4.4)  and  let 

P(A)  =  ABH-(1- A)Pi, 
C(A)=:AC+(1-A)C,-, 


where  Bi  and  Ci  are  matrices  defining  the 
initial  problem  at  A  =  0,  and  correspond  to 
the  par8Lmeter  vector  a  in  Theorem  3.3.  The 
structure  of  the  homotopy  map  /?(A,x,a)  for  the 
input  normal  form  is  now 
Fi{X,X,  a)  =  {PTiBiX)  +  P2Bm)V  +  MeBmV, 

(4.30) 

F2{X,  X,  a)  =  R{CmQ2  ~  C'(A)Qi2)  +  RCmMo, 

(4.31) 

where 


denotes  the  parameter  variables  Bi  G 
Ci  €  R'^”  , 


denotes  the  independent  variables  Bm  and  Cm 
corresponding  to  x  in  Theorem  3.3,  and  A,  P,  C, 
V,  R  are  constants  as  in  Section  2. 

The  Jacobian  matrix  of  p(A,a:,a)  has 
nm^n  -f  riml  rows  and  (um  +  n)(m  +  /)  +  1 
columns.  Since  Pi(A,a:,a)  does  not  involve  C,- 
and  F2{X^x^a)  does  not  involve  P, 

i>c,Pi(A,x,a)  =  0,  PB,P2(A,x,a)  =  0, 
The  Jacobian  matrix  is 
Dp(X,x,a)  = 

/  Da  Pi  Db^Pi  Dc^Pi  Db^Pi  0  \ 

\DxF2  Db^Pz  0  DciF2  J 

(4.32) 

The  following  lemma  will  be  used  in  the  proof 
of  Theorem  4.5. 

Lemma  4.4.  Let  .4  ,  P,  C,  Ajy  Bj^  Cj,  P,  Q,  P, 
P/f  Q/i  Pi 2  and  U  be  defined  as  above.  Then 
Qi  =  Qi,  A  =  Pi,  (4.33) 

Qi2  =  QnU-^,  P12  =  PnU,  (4.34) 
Qj  =  /,  Pj  =  n,  (4.35) 

Q2  =  UU'^,  P2  =  U-^W-K  (4.36) 

In  addition,  Pi2,  Qi2i  P12,  and  O12  have  full 
column  rank. 

Proof.  Equations  (4.20)  and  (4.23)  can  be 
written  in  the  form 

{A  0\fQi  Qi2\ 

\0  AmJ\Q'[2  Q2) 


8 


Y.  Wang,  D.  S.  Bernstein,  and  L.  T.  Watson 


C^RC  -Cf^RCm\_n 
+  \-ClRC  C^RCm  J 

Expanding  these  equations  yields 

AQi  +  QiA^  +  BVB'^  =  0,  (4.37) 

AQi2  +  Q12AI  +  BVBl  =  0,  (4.38) 

AmQ^  +  Q2AI  +  Bm  VBl  =  0,  (4.39) 
A'^Pi  +  PiA  +  CFRC  =  0,  (4.40) 

A^Pu  +  PnA„,  -  C/^RC„  =  0,  (4.41) 

AIP2  +  P^Am  +  d^RCm  =  0.  (4.42) 

Comparing  (2.7a)  with  (4.37),  and  (2.7b)  with 
(4,40)  yields  (4.33). 

If  the  definitions  >1,7,  =  U'^^AmUy 
:=  and  Cm  =  CmU  in  Theorem  4.3 

are  substituted  into  (4.17)  then  (4.17)  becomes 
A.J'U'^  +  UU^Al  +  BmVBl  =  0,  (4.43) 

Alu-'^nu-^  +  U-'^ilU-'^Am  +  =  0. 

(4.44) 

Comparing  (2.7a)  and  (2.7b)  with  (4.43)  2uid 
(  l.-M)  yields  (4.36). 

If  .4,r,  =  U  ^ AmU i  Bm  —  U  ^Bmi  and 
Cm  =  Cmf "  are  substituted  into  (4.38)  and  (4.41) 
and  the  resulting  equations  aire  compared  with 
(2,7a)  and  (2,7b),  then  (4.34)  follows.  Comparing 
(4-17)  and  (4,18)  with  (4.39)  and  (4.42)  yields 
(4.35). 

Finally,  since  and  P2  are  nonsingular,  from 
Section  6  in  (Ge  et  a/.,  1996)  it  follows  that 
Qi^y  and  Pi 2  have  full  column  rank.  Since  U  is 
nonsingular,  from  (4.34)  it  follows  that  Q12  and 
Pj2  also  have  full  rank.  Q.  E.  D. 

Thi^yrcni  i  5  Let  P/  and  Q/  be  defined  as  above. 
Then  Dp{\.r.a}  given  by  (4.32)  has  full  column 
rank  for  0  <  A  <  1,  i.e.,  the  homotopy  map 
{4.30)-(  l  31)  is  transversal  to  zero  for  0  <  A  <  1. 

Proof.  To  prove  Dp{X,z,a)  given  by  (4.32)  has 
full  column  rank,  i.e,, 

rank  (D/>(A, x, a))  =  n^m  -f  riml, 
it  suffices  to  prove  that 

rank(DaPi)  =  rank(DB,Pi)  =  (4.45) 

rank(Dap2)  —  rank(Dc,P2)  =  (4.46) 

Since  V  and  R  are  constant  symmetric  positive 
definite  matrices,  without  loss  of  generality  set 
=  /  in  (4.30)  and  P  =  /  in  (4.31)  .  Using 
Lemma  4.1  to  compute  i)BjPi(A,x,a),  ordering 
Bi  and  Pi  by  columns, 

BB.Fi(A,x.a)  =  £)B.(Pi^2B(A)) 

=  {1  -  X)DB.{Pl2Bi) 

=  (1-A);„,®  (4-47) 


Ordering  Q  and  Fj  by  rows  gives 

BciF2(A,*,a)  =  I>c<(-C(A)(5i2) 

=  (A-l)I>c,(aQi2) 

=  (A-l)/,®(5f2-  (4.48) 

Now  finally,  using  Lemma  4.4,  (4.47),  and  (4.48), 
the  rank  statements  of  (4.45)  Euad  (4.46)  follow. 

Thus  the  homotopy  map  (4.30)-(4.31)  for  the 
input  normal  form  psirametrization  of  {Amj  Bm^ 
Cm)  for  the  model  order  reduction  problem  is 
transversal  to  zero.  Q.  E.  D. 


4.3.  TrsLnsversaJity  of  bomotopies  based  on 
Ly^s  formulation 

In  Ly’s  formulation  (Ly  et  aJ.,  1985),  the 
reduced  order  model  is  represented  with  respect 
to  a  basis  such  that  Am  is  a  2  x  2  block-diagonal 
matrix  (2x2  blocks  with  an  additional  1x1 
block  if  rim  is  odd)  with  2x2  blocks  in  the  form 

C  :)• 

Bm  is  a  full  matrix,  and  Cm  =  ((C„)i  (C,„)2  ••• 
(Cm).' •••  (Cm)r),  where 


(Cm)r  “(1  ^  ^)  >  Um  is  odd. 


Let  S  be  the  set  of  indices  of  those  elements 
of  Am  which  are  independent  variables,  i.e., 
5=  {(2,1),  (2,2),  ...,  (2z-,2f-l),  (2i,20,  ..., 
(ntn,nm)}.  To  minimize  the  cost  function 
J(Am^BmiCm)i  Consider  the  Lagrangian 

L(/lm.Bm,Cm,(3)  =  tT[QR+{AQ+QA^  +  V)  P], 

(4.49) 

where  the  symmetric  matrix  P  is  a  Lagrange 
multiplier,  Q  satisfies  (4.20),  A,  P,  and  V  are 
defined  in  Section  4.2.  Setting  dLjdQ  =  0  gives 
(4.23);  Q  and  P  are  symmetric  positive  definite 
and  can  be  partitioned  as  in  (4.22)  and  (4.24).  A 
straightforward  calculation  shows 


dL 

d{Am)ij 

dL 

dBm 

dL 

d(Cm)ij 


=  2(F^Qi2  +  P2Q2)fj,  {i,j)  €  S, 


=  2{P^:,B  +  P2B,n)V, 

(-Qr2C^BCm) 


+  tr  (Q^ClRCm)] 


=  2R{CmQ2-CQu).j,  i>l. 

(4.50) 

Let 


A{X)  =  A,  B(A)  =  AB  +  (l-A)Bi, 
C(A)  =  AC+(1-A)C<, 
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where  3%  and  Ci  play  the  same  role  as  in  Section 
4.1.  Let 

=  (A^2Qi2  +  P2Q2), 

Hc,^{\x,Ci)  =  2^^  ~  R{CmQ2  -  C'(A)Qi2), 

(4.51) 

where  in  Ha^  only  those  elements  corresponding 
to  the  independent  variables  of  Am  sure  nonzero 
and 

/  (-<4^)5  \ 

x^l  Vec  (Bm)  (4.52) 

\Vec  {Cm)rJ 


denotes  the  independent  variables,  (>1^)5  is  a 
vector  consisting  of  those  elements  in  Am  with 
indices  in  the  set  «S,  i.e., 

{Am)s  =  ((^m)21)  (^m)22)  *  •  *>  (-4m)n«nm)  i 


[Cm)r  is  the  matrix  obtained  from  rows 

T-  =  {2,...,/}  of  Cm. 

The  homotopy  map  p(A,x,a)  for  Ly’s 

formulation  is  now  defined  as 

Fi(\,x,a)=[HAjX,x)]^,  (4.53) 

F2{X,z,a)  =  Vec  [HBjX.x.Bf)],  (4.54) 

FA\,x,a)  =  Vec  [//c«(A,r,C.)]^  ,  (4.55) 

where  again  the  subscripts  S  and  T*  select  the 
appropriate  matrix  elements,  and 


denotes  the  parameter  variables.  As  discussed  in 
Section  4.2,  without  loss  of  generality  set  V  =  / 
in  (4.54)  and  /i  =  /  in  (4.55). 

The  Jacobian  matrix  £)p(A,  x,a)  of  p(A,x,a)  is 

/DxFx  DrFi  0  0  \ 

DxF2  D,F2  Db.F2  0  .  (4.57) 

\DxF3  D^F^  0  Dc^FsJ 


Lemma  4.6.  Suppose  rank(DxFi)  =  rim-  Then 
the  Jacobian  matrix  (4.57)  has  full  column 
rank  for  all  0  <  A  <  1,  i.e.,  the  homotopy 
map  (4.53)-(4.55)  is  transversal  to  zero  for  zJl 
0<  A<  1. 

Proof.  A  similar  proof  to  that  in  Section  4.2 
yields 

rank(D5,F2)  =  rnrim  for  A  1.  (4.58) 


Ordering  C,*  and  F3  by  rows  gives 
L>CiF3(A,^,a)  —  Dci(“C'(A)Qi2)r 
=  (A  —  l)I>Ci(C',Qi2)r. 

=  (A  -  l)Dci[{Ci)r  Qi2] 


1  times 

_ _  - _ 

/o 

Ql2 

0  .. 

•  ®  ^ 

0 

0 

QJ2  •• 

.  0 

\0 

0 

0  .. 

•  Qi2^ 

(4.59) 

and  then  as  before 

rank(2?c,F3)  =  (/-l)nm  for  A^l.  (4.60) 
Note  that 

rank(DxFi)  = 

which  completes  the  proof.  Q.  E.  D. 

Note  that  there  are  only  components  in  Fi 
but  (/  -f  m)nm  +  1  independent  variables  in  x 
and  A.  As  /  +  m  1  usually  in  real  problems 
which  have  been  considered  previously  (Ge  et 
aJ.,  1996),  all  Jacobian  matrices  of  Fi  in  those 
problems  satisfied  the  full  rank  condition.  Since 
each  of  Quy  P12,  Q2,  and  P2  are  implicit 
functions  of  x  and  A(A),  and  one  can  not  give 
explicit  expressions  for  D^Fi  or  Da,Fi  as  in 
(4.59)  for  Dc^Fz  (which  show  clearly  the  rank 
conditions),  it  W2is  necessary  to  assume  that 
rank(DxFi)  =  Um  in  Lemma  4.6.  To  guarantee 
the  full  rank  of  Dp  without  this  assumption, 
instead  of  using  (4.53),  let  x  —  (ry,C),  V  € 

F,(A,x,a)  =  a[h^„(A,i)]^  +  (1  -  A)(»7  -  7/0). 

(4.61) 

with  rim  independent  parameter  variables  in  r/o, 
which  gives 

£>^,Fi  =  (l-A)/„^  for  A^l.  (4.62) 

Combining  (4.58),  (4.60),  and  (4.62)  completes 
the  proof  that  the  map  (4.61),  (4.54),  and 
(4.55)  is  transversal  to  zero.  Note  that  the 
homotopy  construction  in  (4.61)  is  a  theoretical 
convenience,  and  in  practice  the  choice  (4.53)  has 
been  entirely  satisfactory. 
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5,1,  Counterexample  for  optimal  projection 
bomotopies 

The  zero  set  Pa^(O)  of  a  given  homotopy 
map  based  on  the  optimal  projection  equations 

(4.1)-(4.3)  is  not  cdways  bounded,  as  shown  by 
the  following  2-dimensional  example. 

The  system  (Kabamba,  1985b)  is  given  by 
,  /-0.25  -0.4^ 

"^“(^-0.4  -0.72/’  1,1.2)’ 

C=(l  1.2).  (5.1) 

For  the  system  (2.1)-(2.4)  defined  by  (5.1),  the 
solution  set  of  the  optimal  projection  equations 

(4.1) -(4.3)  contains  an  isolated  solution  and  a 
one-dimensional  manifold  of  solutions. 

The  isolated  solution  of  this  system  is 
>1^  =  (-0.838521),  Bm  =  (1.537575), 

Cm  =  (1.537575), 

which  was  obtained  by  both  POLSYS  from 
HOMPACK  (Watson  et  ai.,  1987)  and  by  a 
homotopy  approach  (Zigic  et  aJ.,  1993b).  The 
one-dimensional  manifold  of  solutions  can  be 
derived  directly  from  equations  (4.1)-"(4.3)  as 
follows. 

Let  IV,  =  (^‘),1/,  =  (ui,U2),E  =  (r,  K  =  /, 

and  R  —  I.  The  optimal  projection  equations 

(4.1) ~(4.3)  for  this  problem  can  be  written  as 
0  =  — 0.25u>iUi£r  —  {)Aw\W2Uia‘  —  0,4wiU2(t 

—  0.72wiW2U2(t  —  0.2bwicr  —  0.4i£;2C^ 

-}-  til  +  1.2u2, 

0  =  — 0,25uiUJ2tii<r  —  {)Aw\\x\(t  —  ^Aw\W2U2C 

—  0.72u  —  0.4u;i(7  —  0.72u;2^^ 

+  1.2tii  +  1.44ii2t  (5.2) 

0  =  — 0-25u'itij(T  —  0.4u’2tii<T  —  0.4u;itiiti2<^ 

—  0.72u>2tii ti2^  ~  0.25tiiO’  —  OAU2O’ 

-h  wi  4-  I .2u'2, 

0  =  — 0.25 u;i til U2fr  “  ^Aw2’U\U2(t  —  0.4u;iti2^ 

—  0.72u>2ti2^  —  0.4ui£r  —  0.72ti2<^ 

-h  1.2u;i  -h  1.44ti;2, 

0  =  u^itii  -h  u^2ti2  —  1- 

The  triple  (A  m  1  Bm  I  Cm )  Is  given  by 

=  TAG^  =  (^2) 

=  u;i(— 0.25ui  —  O.4U2) 

+  u;2(-0.4iii  —  0.72ti2)j 

=  TB  =  (u,  U2)  ( 1^2)  = 

C7„  =  CG^  =  (1  1.2)  =«;i  +  1.2tx;2, 

where  F  =  C/j  and  G  =  , 


The  zero  set  of  (5.2)  contains 

{{WuUuE)  :  wi  =  -1.2u;2,  =  -1.2ti2, 

1 


ti2  = 


2.44u;2’ 


(7  =  0} 


which  is  unbounded.  Every  point  in  this  set 
corresponds  to  the  same  triple  (Am,  Bm,  Cm)’ 


Am  =  -.0491803,  Bm  =0,  Cm  =  0, 

The  homotopy  map  based  on  the  optimal 
projection  equations  is 

Ui  >1(A)  Wi  E  +  E  wf  A^{X)  +  Ui  BVB^  =  0, 
A’’ (A)  E  +  E  Ui  >l(A)  Wi  +  C^RC  Wi  =  0, 
UiWi-I  =  0,  (5.4) 


where  A(A)  =  AA  -h  (1  -  A)!),  and  D  is  part  of 
the  parameter  vector  a  in  Theorem  3.3.  The  zero 
set  />J^(0)  of  this  homotopy  map  for  the  system 

(5.1)  includes  the  subset 

{(A,  Wx.UuE):  0  <  A  <  1,  u;i  =  -^l.2w2, 

«1  =  -I.2ii2,  “2  =  2;^,  <r  =  0},  (5.5) 

which  is  unbounded.  This  example  shows  that 
the  zero  set  Pa^(O)  of  a  homotopy  map  can  be 
unbounded  and  yet  some  zero  curves  may  still 
converge  to  isolated  solutions. 

Note  that,  in  pr2u:tice,^^e  algorithm  in  (Zigic 
et  ai.,  1993b)  always  maintains  rank(E)  = 
where  =  1  in  the  above  example.  Solutions 
with  E  =  0  in  the  above  example  never  come 
into  play.  Boundedness  of  ^(0)  for  the  optimal 
projection  equations  (4.1)-(4.3)  can  indeed  be 
guar2Lnteed  with  more  sophisticated  mathematics, 
a  slightly  different  homotopy  map  from  the  one 
used  in  practice,  and  complex  arithmetic  for  the 
curve  tracking.  This  is  pursued  in  Section  5.3. 


5.2,  Simplification  and  example  for  input 
normal  form  homotopy 

The  following  corollary  is  needed  to  simplify 
the  homotopy  map  based  on  the  input  normal 
form  formulation  for  the  optimal  model  order 
reduction  problem. 

Corollary  5.1.  Let  Aj^  Rj^  Vj  he  defined  as  in 
Section  4.2,  partitioned  as  in  (4.21),  let  Am 
be  stable,  and  Qj  satisfy  (4.20).  To  minimize 
(4.19)  under  the  constraints  (4.17)  and  (4.20), 
the  following  two  Lagrangians  are  equivalent: 

Ll{Am,  Bm,  Cm,  Qj,  ^c,  Mo,  Pi ) 

=  tr  +  {A^  "b  PmV Pm) Me 

+  +  ^Am  “h  C^RCm)Mo 

+  {AjQi-\-QiA^i  ■¥Vi)Pi\,  (5.6) 
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where  the  symmetric  matrices  Af©,  Me,  and  Pj 
are  Lagrange  multipliers  introduced  in  Section 
4.2,  and 

L2{A,r„Bm,Cm,Ql,Pl)  =tT[QlRl 

AiAiQi  +  QiAj  ■\-Vi)Pi\,  (5.7) 

where  Qj  is  restricted  to  the  form 


the  Lagrange  multiplier  Pj  is  restricted  to  the 
form 


and  fi  =diag(u;i,  •  •  is  a  positive  definite 

matrix. 


Proof.  The  proof  is  straightforward.  Setting 
()i/dQi  =  0  gives  the  same  equation 

AjP,  +  PiAi  +  Ri  =  0  (5.8) 

111  both  cases.  Expanding  (4.20)  and  (5.8)  yields 
th*  equations  for  and  P2.  In  the  first  case 
ArnQl  +  +  PmV ~  Oj 

AlP2  +  P2A„  +  ClRC„=Q. 

Siiii  *'  the  constraints  (4.17)  and  (4.20)  should  be 
.s.iusfifd  and  Am  is  stable,  it  follows  that  at  a 
constrained  minimum 

O2  =  7n«.  P3  = 

Q.  E.  D. 


liif  [lartial  derivatives 
can  computed  a.s 


dL2 

dBm 


tind 


dLi 

dCm 


dL-2 

tilim 


'2{Pl2B  +  nBm)V, 


of  Lj 


f>/-s 

dCm 


2R{Cm-CQl2). 


'I'he  corri'spf)nding  homotopy  map  (4.30)  and 
(4.31)  i.s  now  simplified  as 


A>(A,x,a) 


(\ec{HB,(X,x,a))\ 

\Vec{Hc^{\,z,a))J' 


where 

//B.(A,x,a)  =  (Pi^,B(A)  +  nS„)V, 
IIc^{X.x,a)  =  R{Cm  -  C(X)Qn). 


The  zero  set  pj'(O)  of  a  homotopy  map  based 
on  the  input  normal  form  formulation  given  by 
(Ge  et  aJ.,  1994)  is  not  always  bounded,  as 
shown  by  the  following  2-dimensional  example. 
The  system  is  given  by 

/-0.895116  0.612237  \  o  ^ 

^  0.612237  -0.447393  j  ’  V  1  y  ’ 

C  =  (-2  1).  (5.9) 


According  to  (Ge  et  al,,  1994),  the.  initial  point 
and  the  triple  (A(A),  B(A),  C(A))  are  chosen  as 
follows: 

1)  Transform  the  given  triple  (A,  B,  C) 
to  bzJanced  form  (^6,  Ct),  such  that 

Ab  =  T“MT,  Bb  =  r“^B,  and  Cb  =  CT  satisfy 

0  =  A^A  +  AA^  -h  BbV BJ  , 

0  =  A^ A  +  AAft  “h  RCb, 

with  a  positive  definite  matrix  A  =  diag  (di,  d2> 

*  *  'i  ^n)j  di  >  di^i* 

The  bsJztnced  form  of  (5.9)  is 

.  -0.25297  -0.5  \  -1.232  \ 

_o.5  -1.0896  y  ’  “  1^-1.866  J  ’ 

C6  =  ( -1.232  -1.866), 

with 


/  0.866 

0.5  A 

A=r 

®  ) 

V  0.5 

-0.866  J  ’ 

'  1  0 

1.5978  J 

2)  For  Tim  —  1)  the  parametrization  (A(A),  P(A), 
C(A))  is  chosen  as 

A{X)  =  Ayl  +  (1  -  A)yl,-  =  ) 

_  /-0.6422A- 0.25297  0.612237A  ^ 

0.612237A  0.6431A-  1.08967  ’ 

P(A)  =  AB  +  (l-A)B,=  (^;{^j) 

_  ^-1.232 -0.768A^ 

C(A)  =  AC  +  (l-A)C,  =  (ci(A)  C2(A)) 

=  (-1.232 -0.768A  A)  =  5"^  (A). 


where 


4,=  ( 

,  _  7-1.232A 

•  ~  I  0  ;  ’ 


-0.25297  0 

0  -1.0896 


C,  =  (-1.232  0). 


For  brevity,  ai(A),  a2(A),  a3(A),  6i(A),  62(A), 
ci(A),  and  C2(A)  will  be  denoted  by  ai,  02,  os) 
61,  62,  Cl,  and  C2  respectively  in  the  following. 
As  discussed  in  Section  4.2,  without  loss  of 
generality,  set  K  =  /  and  R=  1. 

For  any  0  <  A  <  1,  Bm  E  R>,  Bm  ^  0,  let 

—  r2 

Am  =  Cm  =  -VnBm, 

5  r  M(6i- 62)61  12 

^  Lai  +  Am  —  Afa2J 

0261  —  62(61  —  Am) 


M 


{Pn)i2  = 


{Pi2)n  = 


61(03  -f  Am)  -  6202* 

Gm(g26l  —  0162  —  Am62) 


02  —  (oi  +  Am)(u3  +  Am)  ’ 
hCm  “  (^12)12(^3  H-  Am) 


02 
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\  (^12)11  ,/s  X  _  (^’12)12 

“Cl"’ 

Then 

p(\,x,a)  =  0, 

i/(A)Q/  +  Q/i?’(A)  +  F/(A)  =  0, 

Aj{X)Pi  +  P/i/(A)  +  RiiX)  =  0 

are  satisfied.  The  zero  set  pJ^(O)  of  this 
homotopy  map  includes 

{(A,  Bfn ,  Cm)  •  0  <  A  <  1,  Cm  ^ 

(5.10) 

Clearly,  (5.10)  is  unbounded.  If  Bm  0,  then 
Am  is  stable,  [Am,  Bm)  is  controllable,  and  (ylm, 
Cm)  is  observable. 

5.3.  Homogeneous  transformation  to  avoid 
solutions  at  infinity 

As  shown  by  the  examples  in  Sections  5.1 
and  5.2,  the  polynomial  systems  (4.1)-(4.3)  or 
(•1  30) -(4.31)  may  have  solutions  at  infinity,  and 
contains  paths  that  diverge  to  infinity  as 
A  ajiproaches  1.  Solutions  at  infinity  can  be 
avoulf'd  via  the  following  trainsformation  (Morgan 
and  Sommese,  1989),  (Morgan  and  Sommese, 
19S7a).  which  will  be  used  in  Section  5.4. 

Lf't  /(c)  =  0  be  a  polynomial  system  of  N 
equations  in  N  unknowns,  where  z  E  C^,  and 
define  f'{z*)  as  the  homogenization  of /(z): 

fji-')  =  W^o), 

(6.11) 

where  dj  =  deg(/j).  /'(z')  =  0  is  a  system  of  N 
homogeneous  equations  in  *f  1  unknowns. 

Note  that,  if  /'(z°)  =  0,  then  /'(cr°)  =  0 
for  any  complex  scalar  c.  Therefore,  we  may 
take  “solution.s*'  of  /'(c')  =  0  to  be  (complex) 
lines  through  the  origin  in  The  set 

of  these  lines  is  called  complex  projective 
space,  denotf'd  by  a  smooth  compact 

A'-cornplex-dimensional  manifold.  It  is  natural 
to  view  as  a  disjoint  union  of  points 

[(co, . . .  t  2,v)]  with  zo  0  aiid  the  “points  at 
infinity'’,  the  points  [(roi •  •  •  i^^jv)]  with  Zo  =  0. 
The  solutions  of  /'(z^  =  0  in  P^  are  identified 
with  the  solutions  and  solutions  at  infinity  of 
/(z)  =  0  as  follows. 

First,  the  solutions  to  /(z)  =  0  can  be 
identified  with  the  solutions  to  /'(z')  =  0  with 
zq  0.  Explicitly,  if  L  E  P^  is  a  solution  to 
/'(r')  =  0,  and  z'  E  X/,  with  z'  =  (zq,  . . . ,  z;v)  and 
zo  ^  0,  then  z  =  (zi/zo,  Z2/Z0, . . . ,  z^/zq)  is  a 
solution  to  /(z)  =  0.  On  the  other  hand,  if 
zEC^  is  a  solution  to  /(z)  =  0,  then  the  line 
through  z'  =  (l,z)  is  a  solution  to  /'(z')  =  0 
with  zo  =  1  ^  0.  A  “solution  to  /(z)  =  0  at 
infinity”  is  simply  a  solution  to  /^(z')  =  0  (in 
P^)  generated  by  z'  with  zq  =  0. 


Define  a  homotopy  map  (in  P^) 
h(z^  A)  =  (1  ^  A)7y (zO  +  A/'(z0,  (5.12) 

where  is  a  homogeneous  system  of  iV 
polynomials  in  N  +  1  variables,  and  7  is  a 
randomly  chosen  complex  number.  Intuitively, 
let  be  chosen  so  that  its  homogeneous 
structure  matches  that  of  /'.  Precisely,  let 
S  E  P^  be  the  set  of  common  solutions  of 
/'(z')  =  0  and  g'{z*)  =  0.  Then  for  each  s  E  S 
the  following  conditions  must  hold.  For  s  E  S 
let  K  denote  the  full  connected  component  of 
solutions  of  g^{z')  =  0  with  s  E  K. 

If  s  is  a  geometrically  isolated  solution 
of  g\z^)  =  0,  assume  that:  a)  s  is  also  a 
geometrically  isolated  solution  of  /'(z')  =  0,  and 
b)  the  multiplicity  of  s  as  a  solution  of  g\z^)  =  0 
is  less  than  or  equal  to  the  multiplicity  of  s  as  a 
solution  of  /'(z')  =:  0. 

If  s  is  not  a  geometrically  isolated  solution  of 
g*{z^)  =  0,  assume  that:  a)  K  is  contained  in  5, 
b)  AT  is  the  full  solution  component  of  /'(z')  =  0 
containing  s,  c)  K  is  a  smooth  manifold,  and  d) 
at  each  point  z^  E  K  the  rank  of  V^'(z°)  is  the 
codimension  of  K. 

Let  S'  denote  the  solution  set  of  ^'(z')  =  0 
in  P^  —  S.  Under  these  assumptions,  the  basic 
result  is  the  following  theorem. 

Theorem  5.2  (Morgan  and  Sommese,  1987b). 
Assume  the  points  in  S'  are  nonsingular  solutions 
of  y'(z')  =:  0.  For  any  positive  r  and  for  all 
but  a  finite  number  of  angles  0,  if  7  =  re’^, 
then  h~^(0)  Pi  ((P^  —  S)  x  [0, 1))  consists  of 
smooth  paths  and  every  geometrically  isolated 
solution  of  /'(z')  =  0  not  in  S  has  a  path  in 
(P^  •“  S)  X  [0, 1)  converging  to  it. 

Let 

N 

L(z’) 

i=0 

where  6^  ^  0  for  some  i. 

UL^meP^  \L{z'):^0} 

is  the  Euclidean  coordinate  patch  on  P^  defined 
by  L.  Note  that  [//,,  which  is  an  open  dense 
submanifold  of  P^,  can  be  identified  with 
via 

[(zqi  •  •  •  1  z/v^)]  y  “  (zq ,  •  •  •  5  Zj_i,  z,'.j-i, . . . ,  z^), 

L[z  } 

where  6,-  ^  0. 

The  following  theorem  from  (Morgan  and 
Sommese,  1987b)  shows  how  to  keep  the 
homotopy  process  in  complex  Euclidean  spaoe, 
even  though  the  basic  theorem  is  formulated  in 
P^. 
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Theorem  5.3.  Assume  the  points  in  are 
nonsingular  solutions  of  g'{z')  =  0.  Then 

A-i(0)n((P^-5)x[0,  l))  CUlx  [0, 1], 

for  almost  all  Ul  and  all  but  a  finite  number  of 
angles  6. 

For  computations,  the  coordinate  patch  Ul 
is  realized  via  a  projective  transformation  as 
follows.  With  homogeneous  h  in  the  variables  z,* 
for  i  =  0  to  Ny  let 

N 

2o  =  Yl0iZi-\- Pa,  (5.13) 

1  =  1 

where  the  pi  are  constants  and  /?,*  ^  0  for  all  i. 
The  projective  transformation  of  h  is  the  system 
//  of  N  polynomials  in  the  N  variables  z,*  for 
1  =  1  to  where  Hj  =  hj,  with  (5.13)  defining  zq 
in  terms  of  the  other  variables.  By  Theorem  5.3, 
the  homotopy  paths,  including  end  points,  are 
completely  represented  in  via  H.  The  finite 
solutions  of  /(x)  =  0  are  recovered  via  Zj  ^  Zi/zo 
for  1  =  1  to  N.  If  zo  =  0,  then  the  solution 
is  at  infinity.  This  concludes  the  background 
discussion  of  polynomial  system  theory. 

5  i.  fiomogeneous  trajisformation  of  optima] 
projection  bomotopies 

In  this  section  the  homogeneous  transformation 
introduced  in  Section  5.3  is  used  to  prevent 
unbounded  zero  sets  for  optimal  projection 
homo  topics .  Consider  the  polynomial  system 
given  by  (4.1)-(4.3)  sind  the  corresponding 
optimal  projection  homotopies  defined  in  Section 
4.1  The  start  system  at  A  =  0  is  taken  as 

r-,.i(0)ir,En7  +  eic7'>i(0)^  +  CiB(0)  =  0, 

A(0fVl'^  +  Ul'i:UiA(Q)Wi+C{0)Wi  =  0, 

^  iH'i  - =0.  (5.14) 

where  A(0)  =  D  ==  A  —  cln,  €  is  a  constant, 
/1(A)  =  A.4  (1  -  A)D.  The  target  system  (  at 

A  =  1)  is  (4.1)-{4.3). 

According  to  Section  4.3,  the  homogenization 
of  the  target  system  (4.1)-(4.3)  can  be  taken  as 

+  zli:'w(^A'^  +  zIu[bvb'^  =  o, 

zlA^U['^'L'  +  u(^IfU[AW[  +  zlCf^RCW[  =  0, 

U[W[-zllr,^=^Q,  (5.15) 

where  z  =  (vec(f/i),  vec(kFi),  vec(S)), 

[/{(zq,  . . . , zyv)  =  ztiU\(z\/zoy...,zr^/zQ), 
W[(zo,...,zn)  =  zoWi{zilzQ,...,zslzo), 

E  (zq,  .  .  .  ,  Zat)  =  ZoE(zi/zo,  ...»  Zf^jzo). 

The  corresponding  homogenization  of  the  start 
system  is 

U[DW[T.'WC  +  +  zlU[Bi  =  0, 

zlD^U[^T.'  +  ufT.'U[DW[  +  zlCiW[  =  0, 

V[W[-zllr,^  =  Q,  (5.16) 


where  Bi  =  5(0)  and  Ci  =  C(0). 

Theorem  5.4.  If  C,*,  ajad  e  can  be  chosen 
such  that  (5.16)  and  (5.15)  have  no  common 
Zo  0,  E'  ^  0  solutions,  and  all  zq  7^  0,  E'  0 
solutions  of  (5.16)  are  nonsingular,  then  every 
geometrically  isolated  solution  of  (5.15)  has  a 
path  in  converging  to  it. 

Proof.  If  e  =  0,  (5.15)  and  (5.16)  have  the  same 
Zo  =  0  solution  set  (corresponding  to  solutions  of 
(4.1)~(4.3)  at  infinity).  Since  Bi  and  Ci  can 
be  chosen  such  that  (5.15)  and  (5.16)  have  no 
common  zo  7^  0,  E'  ^  0  solutions  and  all  zq  7^  0, 
E^  ^  0  solutions  of  (5.16)  are  nonsingular,  then 
all  the  conditions  of  Theorems  5.2  and  5.3  hold. 
For  each  point  in  5',  the  associated  path  in 
can  be  tracked  from  A  =  0  to  A  =  1.  This 
will  yield  the  full  list  of  geometrically  isolated 
solutions  to  H{zy  1)  =  0.  No  paths  diverge  to 
infinity. 

If  e  7^  0,  5(A)  =  5^5^,  and  C(A)  =  Cf^RC 
for  0  <  A  <  1  as  in  (Zigic  et  ai.,  1993b),  using  the 
fact  U[W{  =  0  (when  zq  =  0),  it  is  clear  that  the 
Zo  =  0  solution  set  of  (5.16)  is  the  same  as  that 
of  (5.15).  Similarly,  (5.15)  and  (5.16)  have  the 
same  zq  ^  0  solutions  when  E'  =  0.  Note  that 
this  case  corresponds  to  the  counterexample  of 
Section  5.1.  Tcike  S  be  all  the  zq  =  0  solutions 
and  any  solutions  corresponding  to  zq  7^  0  and 
E'  =  0.  Now  e  can  be  chosen  such  that  (5.15) 
and  (5.16)  have  no  other  common  solutions 
and  all  other  zq  9^  0  solutions  of  (5.16)  are 
nonsingular.  Then  the  technical  assumptions  of 
Theorem  5.2  can  clearly  be  met  for  the  common 
solution  set  S.  Thus  Theorem  5.2  and  Theorem 
5.3  hold  for  the  start  system  (5.16)  in  this  case 
(f  9^  0)  also.  Q.  E.  D. 

The  import  of  this  result  is  that  the  real 
solutions  of  (4.1)-(4.3),  which  satisfy  the  rank 
condition 

rank(Wi)  =  rank([/i)  =  rank(E)  =  nm, 
if  they  exist,  must  be  connected  to  the 
solutions  of  (5.16)  in  P*^  —  S.  Technically,  this 
is  guaranteed  only  with  a  complex  multiplier 
7  in  (5.16),  and  only  if  complex  arithmetic  is 
used  and  the  homotopy  curve  tracking  is  done  in 
P^.  However,  all  this  has  never  been  necessary 
in  practice  (Zigic  et  aJ.,  1993b).  Furthermore, 
observe  that  the  solution  set  (5.15)  includes  all 
solutions  with  rank  E'  <  and  thus  one  is 
guaremteed  of  finding  a  reduced  order  model  of 
order  no  greater  than  Um-  Since  (5.15)  represents 
the  optimal  projection  equations  (4.1)-(4.3)  for 
some  stable  A  (A)  for  every  A,  0  <  A  <  1,  it  is 
clear  why  real  arithmetic  suffices  generically. 
Generically,  the  real  solutions  are  isolated,  have 
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constant  rank,  and  vary  smoothly  with  respect  to 
A  (Morgan  and  Sommese,  1989). 

Finally,  for  the  target  system  (5.15),  it  is 
always  possible  to  take  the  starting  homogeneous 
system  as 

Pj4-Qj4  =  0,  j  =  l,...,N,  (5.17) 

where  pj  and  qj  are  positive  constants  such  that 
(5.17)  has  no  common  solution  with  (5.15). 
Since  all  solutions  to  (5.17)  are  nonsingular,  all 
conditions  of  Theorem  5.2  eind  Theorem  5.3  are 
satisfied.  The  drawback  is  that  the  starting 
system  (5.17)  is  totaJly  unrelated  to  (5.15), 
requires  complex  arithmetic,  and  may  take  more 
steps  to  converge. 

6.  CONCLUSIONS 

Probability-one  homotopy  methods  were 
considered  for  the  problem  of  model 

reduction.  The  crucial  requirement  of 
transversality  was  verified  for  sever2J  homotopy 
maps  including  the  pseudogramian  formulation 
of  the  optimal  projection  equations  as  well  as 
variations  based  upon  canonical  forms.  These 
results  guarantee  good  numerical  properties  in  the 
computational  implementation  of  probability-one 
homotopy  algorithms.  Counterexamples  to  the 
hounded  ne.ss  requirement  of  probability -one 
homotopy  theory  were  provided  for  the 
pseudogramian  formulation  of  the  optimal 
projection  equations  and  for  some  formulations 
based  upon  canonical  forms.  Since  a  solution 
may  not  exist  in  any  particular  canonical  form, 
these  results  are  sharp  for  canonicaJ  forms,  where 
unbounded  ness  corresponds  to  nonexistence  of 
solutions.  However,  for  a  reformulation  of  the 
pseudogramian  optimal  projection  equations  in 
complex  projective  space  using  homogeneous 
transformations,  the  boundedness  assumption 
holds  and  thus  global  convergence  of  the 
homotopy  algorithm  to  a  solution  (in  complex 
projective  space)  is  guaranteed.  Both  the 
genericity  of  real  solutions  and  considerable 
computational  experience  (Zigic  et  a/.,  1993b) 
indicate  that  real-valued  homotopies  are  effective 
in  practice  and  thus  it  is  not  necessary  to  track 
the  homotopy  zero  curves  in  complex  projective 
space. 
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Chapter  1 

Introduction 


The  Robust  Fixed-Structure  Control  Toolbox  is  an  integrated  collection  of 
MaTLAB  functions  that  can  be  used  to  synthesize  fixed-structure  controllers 
that  are  optimal  with  respect  to  a  given  performance  measure  and  at  the  same 
lime  satisfy  stability  and  robustness  constraints.  The  Robust  Fixed- Structure 
Contml  Toolbox  is  designed  to  solve  a  large  class  of  problems,  including  de¬ 
centralized  compensation,  reduced  order  compensation,  controller  design  for 
multiple  plant  configurations,  and  real  parameter  model  uncertainty.  The  flex- 
il>ility  of  the  toolbox  routines  is  due  to  the  use  of  a  decentralized  static  output 
feedback  framework  (Chapter  2)  in  problem  formulation,  which  encompasses  all 
of  the  above  problems,  and  allows  them  to  be  solved  with  a  common  solution 
algorithm. 

Oiu'e  a  control  synthesis  problem  ha>  been  transformed  into  the  decentralized 
output  feedback  framework,  the  next  problem  is  to  optimize  the  free 
paramet<Ts  in  the  controller  witli  n‘>p»-et  tr>  one  of  several  performance  criterion 
that  are  included  in  the  Robtist  FtJtd-Strttrtun  ('onfrol  Toolbox  (Chapter  3). 
A  m«>difi#a!  fjuasi-Newton  unconstrained  optimization  algorithm  [4]  is  used  to 
ac.  ( uHpli^h  this. 

I  allowing  the  di.scussion  on  performance  (Titerion  are  demonstrations  of  the 
I raiC'format ion  of  standard  fixed-structure  control  synthesis  problems  into  the 
<l*r»*ntralized  static  output  feedback  framework,  along  with  Matlab  sessions 
showing  how  these  problem>  ar**  srt  up  and  solved  using  the  Robust  Fixed- 
Strvetutr  Control  Toolbox  (Chapter  4). 

A  discussion  of  the  homotopy  algorithms  utilized  in  this  toolbox  can  be 
found  in  Chapter  5,  followed  by  demonstrations  of  formulations  for  the  model 
order  reduction  and  controller  syntlu'sis  problems  (Chapter  6). 

Descriptions  and  syntax  for  all  core  commands  in  the  Robust  Fixed- Structure 
(\jntrol  Toolbox  are  provided  in  this  <io<'ument  (Chapter  7).  Also  included  are 
several  supplementary  function,s.  some  related  to  automating  the  process  of 
transforming  standard  synthesis  problems  into  the  decentralized  static  output 
feedback  framew'ork,  and  others  which  can  be  used  for  generating  initial  con¬ 
trollers  of  a  given  structure  to  optimize.  The.se  supplementary  routines,  while 


1 


2 

not  exhaustive  or  complete,  can  be  used  to  solve  several  commonly  occurring 
synthesis  problems. 


Chapter  2 

Decentralized  Static  Output 
Feedback 


Tliis  section  reviews  the  decentralized  static  output  feedback  problem  for¬ 
mulation  [8.  3]  for  fixed-structure  controller  synthesis.  For  both  continuous  and 
discrete-time  problems,  consider  the  (;n 1)- vector-input,  (/n+p-f  1)- vector- 
out  put  decentralized  system  shown  in  Figure  2.1,  and  define 


y\ 

,  y  = 

ym 

■  di  ■ 

^1 

d  = 

,  6  = 

.  . 

_  . 

tin*  system  G  have  the  realization 


■  A 

B..  1  Bj  1  Bu  1 

Cy 

Pyu\Pyd\Vyu 

C, 

^  C: 

1  P.-u\P.d\P:u  . 

(2.1) 


(2.2) 


2.1  Continuous-Time  Decentralized  Static  Out¬ 
put  Feedback 

In  a  continous-time  framework,  the  realization  of  G  (2.2)  represents  the  dynam- 


i{t)  =  Ar{t)-^Puu(t)  + Bdd{t)  +  Bu.w(^)^  (2-3) 

with  measurements 

y{t)  =  C^X(l)-i-VyuU(t)  +  Vyjd{t)+VyuMi), 


3 


(2.4) 
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Figure  2.1:  Decentralized  Static  Output  Feedback  Framework 
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model  error  outputs 

e(t)  =  CeX{t)  + 'Peu^i^)  A  Dedd{t) +'De^w{t), 
and  performance  variables 

z{t)  =  C,xit)  +  V,uuii)+V,dd{t)+V,^w{t). 


(2.5) 


(2.6) 


Plant  uncertainty  is  represented  as  decentralized  static  output  feedback  by 
the  relations 

=  (2.7) 

where  the  uncertain  matrices  A'  are  not  necessarily  distinct.  To  represent  de¬ 
centralized  static  output  feedback  control  with  possibly  repeated  parameters, 
we  consider 

«.(0  =  ^';.v.(0.  (2.8) 

where  the  matrices  are  not  necessarily  distinct.  Reordering  the  variables  in 
(2.7)  and  (2.8)  if  necessary,  we  can  rewrite  (2.7),  (2.8)  as 


d{l)  = 

ii{l}  =  Ky{t). 


where  A  and  K  have  the  form 


A  =  bio 'k-diag  (/,,;.  A, . 4,  ©A,), 


A."  =  block-dias  ( ,  k  ] 


.  lo..  V  K,.)  . 


(2.9) 

(2.10) 


(2.11) 

(2.12) 


T*-  7  tli<*  number  of  distnu  t  A,  t  is  the  number  of 

r»  t  It  It  Hi'' of  uncertainty  A,,  i  tie  injinb*r<»l  tj;ainsA.,  and 

■:  i--  tli»  number  of  repetitn»n>  of  K  ,  Noir  \  hat  Aj . A,y  and  A.  i , .  . . ,  A.i, 

.ir«  ii  't  n^'iv^^sarily  square  matre  e'x  an-1  tliat 


t  ~  \ 

.\1{»  rn  it  i\ej\ .  K  can  be  wnti»  n  ?t> 


N  O,  —  m 


A.'  = 


(2.13) 
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where  Qhij  QRij  are  defined  as 


On  01  xr, 

Oci0i  xc, 

Or202Xr, 

Oc202XC, 

Or,_i0,„i  xr, 

0c,_  1  1  XCi 

0r.(i-l)xr. 

/r. 

Qru  = 

0c,(j-l)xc, 

/c. 

Or,(0,-j)xr, 

Ocj(0,'-i)xc, 

Or, 4.10, +1  xr, 

Oc, 4-10, 4-1  xc, 

Ori.0yXr, 

Oci,0i,  xc, 

For  convenience,  define 

Lk  =  I  -  DyuK, 

and  assume  that  Lk.  is  nonsingiilar.  Furthermore,  assume  that 


(2.14) 


(2.15) 


Vr,i  = 

0. 

(2.16) 

'^yd^'^eu  ~ 

0  for  all  A  of  the  form  (2.11), 

(2.17) 

VexXLJ:^Vyd  - 

0  for  all  K'  of  the  form  (2.12). 

(2.18) 

In  this  case,  the  closed-loop  dyanmies  ar<‘ 

Hi) 

=  (.4 + 

)  j-(/)  +  [l\  +  /),/A£)e„  )  n’(0. 

(2.19) 

--(0 

T 

II 

;)  rlM  ^  (/),„  +  /X-rfAD,.  )  uit). 

(2.20) 

\vh»Tf 

i  i 

A  +  B„KLl\\. 

/<„  =  E„  +E.K-I.l'Dy,r. 

+  P.-„  A.  7.,^  *(  V 

/>,„  ^  P.-,,  +P.-.,A.7.;^'Pyu, 

(2.21) 

iii  = 

l\i  + Ruf^'iVP.i 

/),„  i  n,„  +P,,A'Aj^;'Py.,., 

L ,  +  P,  U  „ . 

IK  i  —  P.-ii  +  Prti^  7, jj^'pyd- 

The  closed-loop  transfer  function 

(•:n  therefore  has  the  realization 

.4  +  ii.iA(\ 

Bxv  +  Bd^Dcw 

(\  + 

Bzxi  +  Dzd^Dew 

(2.22) 


1  he  nominal  closed-loof)  transfer  fimeHon  (s)  is  obtained  by  letting  A  —  0, 
so  that 


CHAPTER  2.  DECENTRALIZED  STATIC  OUTPUT  FEEDBACK 


7 


2.2  Discrete-Time  Decentralized  Static  Output 

Feedback 

In  a  discrete-time  framework,  the  realization  of  G  (2.2)  represents  the  dynamics 
x(i-l-l)  =  Ax{k)  +  Buu{k)  +  Bdd{k)  +  Bu,‘w{k),  (2.24) 

with  measurements 

y(k)  =  Cyx{k)  +  Vyuti{k)  +  Vydd{k)+Vyu,w{k),  (2.25) 

model  error  outputs 

e{k)  =  C,x{f^)  +  Veuu(k)AVedd{k)+Ve^uw{k),  (2.26) 

and  performance  variables 

z(k)  =  C^x{k)  +  V:yu(k)+V,dd(k)+V,u,w{k)-  (2.27) 

Plant  uncertainty  is  represented  as  decentralized  static  output  feedback  by  the 
relations 

d,(A-)  =  (2.28) 

where  again  the  uncertain  matrices  A'  are  not  necessarily  distinct.  To  represent 
decentralized  static  output  feedback  control  with  possibly  repeated  parameters, 
we  consider 

«»(A)  =  A.':y.(A),  . m,  (2.29) 

where  the  matrices  are  not  necessarily  distinct.  Reordering  the  variables  in 
(2.28)  and  (2.29)  if  necessary,  we  can  rewrite  (2.28),  (2.29)  as 

ci(h)  =:  Sf{k),  (2.30) 

u{k)  =  Ak‘y(A).  (2.31) 

uh»rr  A  and  K  have  the  forms  (2  Ih  and  (2  12),  respectively.  Thus,  the 
.ill* mail ve  representation  of  K  (2.13)  stdl  iin|d>  [  sing  t he  same  definition  of  L^: 
i2  1 3 1 .  an<l  making  the  samf‘  assumptniii"  a>  in  tlir  cont inuoiis-time  framework 
1 2  ITo  (2  18).  the  closed-loop  dynamics  an* 

xik  +  1)  =  (.4  -f  xU-f  -t  (/<„  +  RdSlXu)  uik).  (2..32) 

z{k)  =  (f..  -f  ,)  x{k}  ■+  4-  D;rfADfu.)  u'(k),  (2.33) 

\vti<‘r»‘  wo  make  use  of  the  (h'fuiitiou.s  (2.21).  The  discrete-time  closed-loop 
tr.aiisfer  function  Gzb,.,a(s)  therefore  ha.s  the  realization 

G 

and  the  nominal  closed-loop  transfer  function  (7;u  (c)  has  the  realization 


A  +  UdS(\ 
c-  +  b,dM\ 


Bu,  -b  Bd^Dew 
bzu<  +  bzd^be 


(2.34) 


6-,„  (--)  ~ 


(2.35) 


Chapter  3 

Performance  Criterion 


Tlie  Robust  Fixed-Structure  Contwl  Toolbox  ca.n  be  used  to  synthesize  con¬ 
trollers  that  are  optimal  with  respect  to  a  user-chosen  performance  criterion. 
A  performance  criterion  consists  of  a  cost  function  and  possibly  one  or  more 
constraints.  The  cost  function  represents  some  characteristic  of  the  controlled 
>ystem,  while  the  constraints  represent  properties  that  any  feasible  solution  to 
the  Optimization  problem  must  have.  An  example  of  a  cost  function  is  a  norm 
on  a  closed-loop  transfer  function,  while  examples  of  constraints  include  asymp¬ 
totic  stability  of  the  nominal  closed-loop  system  or  robust  stability  with  respect 
\(}  uncertainties  of  a  certain  size  and  structure. 

riie  performance  criterion  options  available  in  the  Robust  Fixed- Structure 
(\nitrol  TooIIk)x  Kve  described  in  the  following  sections. 


3.1  Continuous-Time  7^ ..-optimal  Performance 


(  out inuous-t ime  'W^-^^ptimal  [)erforman<'e  is  obtained  by  optimizing  the  Tin 
torin  tlic  nominal  closed-loop  system(u  tf)  c)  with  respect  to  the  free  pa- 
ram»  trrs  ofthr  controlh'r.  Since  Hi»-<>pittnal  |)erformance  does  not  account  for 
nil'  *  rtamty  in  the  model,  poor  or  even  <lestabilizing  designs  can  result  if  the 
tru'  plant  diffiTs  from  the  nominal  plant. 

If  (s)  is  an  asymploti<‘all\  stable,  strictly  proper  continuous-time  trans- 
h  r  function  with  the  realization 


C:u(^s) 


A 

Bu.  ■ 

C, 

0 

(3.1) 


then  it  can  be  shown  that  |1  (7,„  (s)  lb  <*an  be  expressed  as 

||C,„(..)ll::=trPC,  (3.2) 

wIktc  P  is  the  solution  of  the  rontiniious-time  matrix  Lyapunov  equation 

{)  =  A'^P+ PA  + R.  (3.3) 
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with  R  —  CjCz  and  V  = 

The  H2  norm  is  only  defined  for  stable,  strictly  proper  transfer  functions. 
Thus,  when  using  the  '?f2-optimal  performance  criterion,  the  closed-loop  system 
will  be  constrained  to  be  asymptotically  stable,  and  the  closed-loop  feed-through 
matrix  Dzw  must  be  identically  zero. 


3.2  Discrete- Time  Performance 


Discrete-Time  7^ 2-optimal  performance  is  obtained  by  optimizing  the  discrete¬ 
time  TZo  norm  of  the  nominal  closed-loop  system  {w  to  z)  with  respect  to  the 
free  parameters  of  the  controller  gain  matrices.  As  in  the  continuous-time  case, 
discrete-time  T^o-optimal  performance  does  not  account  for  uncertainty  in  the 
model. 


If  G-u.(c)  is  a  proper,  asymptotically  stable  discrete-time  transfe  r  function 


with  the  realization 


■  A 

Bu.  1 

C-- 

0 

(3.4) 


tlien  it  can  be  shown  that  1|  G{z)  H-j  can  be  expressed  as 


\\l=  tr  PV  A  DJ^.Dzu,, 


(3.5) 


>  the  solution  of  the  discrete-time  matrix  Lyapunov  equ  ation 

P  =  +  R,  (3.6) 


th  />  =  rJG  and  V  = 

[hr  discrete-time  Tin  norm  is  only  chdined  for  asymptotically  stable,  proper 
tran>hT  functions.  Thus,  the  closed- loop  s\>tem  is  constrained  to  be  asymptot- 
n  allv  stable  when  using  the  discrett-t ime  W:»-oi)timal  performance  criterion. 


3.3  Scaled  Popov  Performance 

S<;ah  d  IN>pov  p(‘rforman(*e  js  characteri/ei.!  l>y  robustness  to  uncertainty  in  the 
f)lani  dynamics,  input  behavior,  output  behavior,  or  any  combination  thereof, 
which  can  be  modeled  by  real  but  uncertain  matrices  inserted  into  the  open-loop 
plant 

x(t)  =  (A  -h  A/^A^  jr(<)  +  {Ru  "f  A/z?AbAb)u(0  -h  Bu,w{i), 

y(i)  =  (Cy  +  A/c  Ac  AV’ (0  +  {Ryu  +  A/d  A/)  +  Dyww(t), 

z(t]  =  C,x{t)^D,uU(t)^D:uW{t).  (3.7) 

lb  synthesize  controllers  whicii  will  exhibit  robustness  to  these  perturbations  of 
the  nominal  plant,  define  A  a.s 

A  =  block-diag(A^ ,  Ab,  Ac,  AD). 


(3.8) 
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It  is  assumed  that  A  is  an  element  of  the  set  of  norm-bounded  matrices  , 
defined  by 

=  {a  G  a  :  crjnax(A)  <  7  (^*9) 

where  A  is  the  set  of  matrices  with  the  specified  internal  structure  of  A.  Chap¬ 
ter  4  gives  examples  of  how  to  transform  (3.7)  into  the  equivalent  decentralized 
static  output  feedback  problem  given  by  (2.2) 

Scaled  Popov  performance  produces  controllers  which  guarantee  closed-loop 
asymptotic  stability  for  all  A  in  the  set  A-y,  while  at  the  same  time  minimizing 
a  bound  on  the  worst-case  7^2  norm  of  the  closed-loop  transfer  function  Gzw{j^) 
subject  to  the  uncertainty  A.  Thus,  the  controllers  not  only  have  a  priori  regions 
of  guaranteed  stability,  but  also  have  a  priori  bounds  on  the  T/o  norm  of  the 
closed-loop  system  for  all  perturbations  within  the  given  class  of  perturbations 

The  worst-case  7^2  norm  for  the  closed-loop  system  from  w  to  subject  to 
the  perturbation  A  is  given  by 

sup  II  l|o  =  sup  tr  P^Bu,B'^,.  (3.10) 

wIhto  is  the  unique,  nonnegative  definite  solution  to  the  matrix  Lyapunov 
r<|uation 

0  =  (.4  +  BdACefl\  +  /^a(.4  +  BdACe)  +  (3.11) 

Sparks  and  Bernstein  [10]  have  shown  that  if 

a>\ rn|>t()tirally  stable,  and  th(T('  exists  a  f)ositive-definite  matrix  P  and  two 
niatrn  es  /  >  0,  the  stabthty  r  matrix,  and  \\  ,  the  sca/zng  matrix  which 

(niiimute  with  the  uncertainty  A.  such  that 

r  =  7^  -  \y(\lli  -  >0.  (3.12) 

and 

0  =  Ajr  +  /Viu  +  +  r.  (3.13) 

,v  i  is] r  +  Z(\  +  irCe.4o,  (3  i4) 

then  the  closed-loop  system  is  asymptotically  stable  for  all  A  €  A.,,  and  the 
worst-case  'H2  norm  obeys  the  inequality 

sup  II  G,„.,a(.s)  II?  <  u{P  +  2y-^CjWCe)BB'^.  (3.15) 


Chapter  4 

Formulation 

(Unconstrained 

Optimization) 


Ill  this  section  we  introciiice  s^'vcral  controller  synthesis  problems,  and  pro- 
\  idc  the  equivalent  decentralized  static  output  feedback  problem.  Unless  specif¬ 


ically  stated  otherwise  within  a  section,  a  plant  of  the  form 

i  =  Ax  ^  fhi^  Diu\  (4.1) 

y  =r  f  >  -h  /'n  -e  D-jU  ,  (4.2) 

c  =  I  xX^  I.'r u  Ei) ir.  (4.3) 

f  r  c< .ntinou>-limf’  problem>  or 

x{k‘ -^  \)  -  \x[k  )  Htiik)  4-  Ditt  ik).  (4.4) 

p\k)  —  ( ’xi  A- 1  i  A  I -f  I)-jtt  {k).  (4.5) 

c(A')  /.\  X!  A' I  -^  / -j  u!  A)  4-  /:nir(A*).  (4.6) 

t  r  r*  o  iinu-.  prol)h‘ins  is  c(»nMd»  r«‘d 


4.1  Centralized  Proper  Dynamic  Compensation 

First  consider  a  full-  or  r^^tluc^^d-order  proper,  centralized  dynamic  compen¬ 
sator 

r,.  =  ArJr.-  +  Dcy.  (4-7) 

./  =  ( -.  r  .  +  D,y.  (4.8) 

or.  in  dis<Tete-time 

x.(A-^  1)  =  A.-Xr(k)^  Bcy[k),  (4.9) 

u(k  )  =  6Vxc(A)  4- Dci/(A).  (4.10) 
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Letting  fC  denote  the  partitioned  matrix 


v4c  Be 
Cc  Dc 


(4.11) 


Lk  is  given  by 


Lk  = 


I  ,  0 
-FCc  I-FDc 


Assuming  the  matrix  I  —  F Dc  is  nonsingular,  it  follows  that  L/c  is  nonsingular, 
and  thus  the  closed-loop  system  consisting  of  (4.1)-(4.3),  (4.7),  and  (4.8)  (or 
(4.4)-(4.6),  (4.9),  and  (4.10)  in  discrete-time)  can  be  written  as  decentralized 
static  output  feedback  with  ni  =  »■  =  di  =  1  and  G  given  by 


A 

0 

0 

5|0|Z?i  ■ 

0 

0 

7 

0  |0|  0 

0 

I 

0 

0  |0|  0 

(' 

0 

0 

F\0\D2 

0 

0 

0 

0  10|  0 

0 

0 

E2I0IEO  _ 

.*\s  an  example,  let  us  use  the  Rt/hust  I iTcd^Structurc  Cofitvol  Toolbox  to 
(ie.sign  an  ?f2-optinial  compensator  with  this  structure.  First,  we  enter  the 
standard  control  problem  plant. 

»  A  •  [zeros(5.1).€ye(^);-l  -2  -24  -12  -24  -4]; 

>>  B  *  [zeros  (5  *  1 ) ; l] ; 

»  C  =  [100000]; 

>>  D  »  0; 

>>  Dl  “  [B .zeros (6 , 1 ) ] ; 

>>  D2  «  [0.1] ; 

>>  El  ■  [C:zeros( 1 .6)] : 

>>  E2  -  [0;  1] ; 

>>  EO  ■  zeros (2,2) ; 

wr  griK-rate  ail  iniiiai  l>-r  ih*  optimal  compensator  by  using  a  bal- 

tiwcfxl  truncation  of  the  full-orii^  r  Hro\)Xuu!\\  compensator,  via  rolqg.  We  also 
<Teate  the  matrices  QLll  and  QRU.  acrording  to  (2.14),  which  define  the  struc¬ 
ture  of  K . 


»  nc  =  4; 

>>  np  *  size(A,  1) ; 

»  [Ac,Bc,Cc,Dc,test]  ■  rolqg( A.B.Dl ,C.E1 ,D,D2,E2,E0 ,nc) ; 
>>  kl  *  [Ac,Bc;Cc ,Dc] 


kl  « 
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-0.0157 

-0.3137 

0.0046 

-0.0016 

-0.0763 

0.3137 

-0.3487 

0.0275 

-0.0093 

0.2845 

-0.0046 

0.0275 

-0.2634 

1.0135 

-0.0112 

-0.0016 

0.0093 

-1.0135 

-0.0360 

-0.0038 

-0.0763 

-0.2845 

0.0112 

-0.0038 

0 

»  QLll  =  eye(5) ; 

»  QRll  =  eye  (5) ; 

»  save  init  kl  QLll  QRll 

VVe  now  transform  the  standard  control  problem  into  the  equivalent  decentral¬ 
ized  static  output  feedback  framework  problem  via  the  realization  (4.12). 

»  A  =  [A,zeros(np,nc) ;zeros(nc,np) ,zeros(nc,nc)] ; 

»  Bu  =  [zeros(np,nc) »B;eye(nc) ,zeros(nc,l)] ; 

»  Bw  =  [Dl;zeros(nc,2)] ; 

»  Cy  =  tzeros(nc,np) ,eye(nc) ;C.zeros(l ,nc)]  ; 

>>  Cz  =  [El .Z€ros(2,nc)] ; 

»  Dyu  =  [zeros (nc.nc) ,zeros(nc, 1) ;zeros(l ,nc) ,D] ; 

>>  Dyw  =  [zeros(nc,2) ;D2] ; 

>>  Dzu  =  [zeros(2,nc) ,E2] ; 

»  Dzw  =  [EO] ; 

»  clear  B  C  D  D1  D2  El  E2  EO  nc  np 
>>  who 

Your  variables  are: 

A  Bw  Cz  Dyw  Dzw 

Bu  Cy  Dyu  Dzu 

W*  iif>u  rnt<T  tiif'  remaining  variabl'>  that  an  nfunled  by  optgain  (see  Chap- 
n  r  7  I,  and  savr  them  to  a  file 

>>  kindex  *  1 ; 

>>  kindexkc  »  5; 

»  clear  kindexkc 
»  indexkc  *  5; 

»  indexkr  «  5; 

»  ctype  “  1 ; 

»  N  =  1; 

>>  who 

Your  variables  are: 


A 


Cy 


Dyw 


N 


indexkr 
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Bu  Cz  Dzu  ctype  k  index 

Bw  Dyu  Dzw  indexkc 

»  save  datafile 

We  load  the  file  init.mat  containing  our  initial  value  for  the  controller  and 
datafile  .mat  containing  the  decentralized  static  output  feedback  problem  data, 
and  save  the  combined  data  to  a  single  file,  runfile.mat. 

>>  clear 

>>  load  datafile 
»  load  init 
>>  who 

Your  variables  are: 


A 

Cy 

Dyu 

N 

ctype 

kl 

Bu 

Cz 

Dzu 

QLll 

indexkc 

k index 

Bu 

Dyu 

Dzw 

QRll 

indexkr 

>>  save  runfile 

\Vf  now  clear  the  workspace  and  execute  opt  gain, 

>>  clear 

»  [fmin,info,noits]  »  optgain( 'runfile 'outfile O 
USING  LINE  SEARCH  ALGORITHM 

ITERATE  FUNCTION  VALUE 


0 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


0,2821656l52736E+00 
0.2815935488427E*^00 
0.2685207335632E4‘00 
0.2660287632697E+00 
0.2646826543299E+00 
0.2635998976307E+00 
0. 26332447 15048E+00 
0.2632156677085E+00 
0.2632052791 173E+00 
0.2631913346650E+00 
0.2631636997435E+00 
0.2631177717331E+00 
0. 263084 1561395E+00 
0 . 2630724598839E+00 
0.2630685095928E+00 


14 
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TOTAL  ITERATIONS  =  14 

UNCMND  WARNING  —  INFO  =  1:  PROBABLY  CONVERGED,  GRADIENT  SMALL 
fmin  = 

0.2631 
info  = 

1 

noits  == 

14 

We  load  the  file  outf  ile.mat  and  examine  the  optimized  controller  parameters. 

»  clear 
»  load  outfile 
>>  who 

Your  variables  are: 
fmin  info  koptl 

>>  koptl 
koptl  = 


0.0032 

-0.2990 

0.0067 

-0.0019 

-0.1655 

0.2990 

-0.3610 

0.0282 

-0.0109 

0.1538 

-0.0067 

0.0282 

-0.2634 

1.0135 

-0.0099 

-0,0019 

0.0109 

-1.0135 

-0.0360 

0.0028 

-0.1655 

-0.1538 

0.0099 

0.0028 

-0.2535 

>> 


4.2  Centralized  Strictly  Proper  Dynamic  Com¬ 
pensation 

( ’nn.sid^'r  a  full-  or  reduced-ord'T  stru  tly  proper,  centralized  dynamic  com- 
p<-ii^ator  having  the  realization 

Xr  =  ArXr+DcV,  (4-13) 

u  =  c;.r,.  (4.14) 

Lotting  A.'  denote  the  block-diagonal  matrix 

K  =  block-diag(.4,,  Sc.Cc),  (4.15) 

it  can  be  verified  that  La.-  is  nonsingular.  Hence,  the  closed-loop  system  con¬ 
sisting  of  (4.1)-(4.3)  and  (1.13)  (1.14)  can  be  written  as  decentralized  static 
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output  feedback  with  m  =  t;  =  3,  =  02  —  ^3  —  I7  given  by 


A 

0 

0 

0 

B\0\Di  ■ 

0 

0 

I 

I 

0  jol  0 

0. 

I 

0 

0 

0  |0|  0 

c 

0 

0 

0 

F\0\D2 

0 

I 

0 

0 

0  |0|  0 

0 

0 

0 

0 

0  |0|  0 

0 

0 

0 

£'2|0|£'o  . 

(4.16) 


We  now  demonstrate  how  to  use  the  Robust  Fixed’ Structure  Control  Toolbox 
to  design  a  continuous-time  compensator  with  this  structure.  First, 

we  enter  the  standard  control  problem  plant, 


>>  A  =  [zeros(5, 1)  ,eye(5)  ;-l  -2  -*24  -12  -24  -4]; 

>>  B  =  [zeros(5, 1) ;  1]  ; 

»  C  ~  [1  0  0  0  0  0]; 

»  D  =  0; 

>>  D1  =  [B,zeros(6, 1)] ; 

»  D2  =  [0.1]  ; 

>>  El  =  [C;zeros  (1 ,6)]  ; 

»  E2  =  [0;1]; 

>>  EO  =  zeros(2.2) ; 


Ne\!.  wo  generate  an  initial  guess  fi)r  the  optimal  compensator  by  using  a  bai 
triinration  of  the  full-order  ^-..-optimal  compensator,  via  rolqg.  VVe  also 
.  r.  air  th.’  matrices  QLij  and  QR/./-  a'  conimg  to  (2.M),  which  define  the  .structure 
of  A,'. 


>>  nc  »  1; 

>>  np  *  size ( A ,  1 ) ; 

>>  [Ac ,Bc .Cc .Dc .test]  =  rolqg<A.B,Dl .C,El,D.D2,E2.E0,nc) ; 
>>  kl  *  Ac; 

>>  k2  *  Be; 

»  k3  «  Cc; 

»  Dc 
Dc  “ 

0 

>>  QLll  =  [eye(nc) ;zeros(nc.nc) ;zeros(l ,nc)] ; 

»  QRll  -  [eye (nc) .zeros (nc.l) .zeros(nc,nc)] ; 

>>  QL21  =  [zeros (nc. nc) ;eye(nc) ;zeros(l ,nc)] ; 

>>  QR21  =  [zerosd  ,nc)  .  1  .zerosd  ,nc)]  ; 
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»  QL31  =  [zeros (nc , 1) :zeros(nc, 1) : 1] : 

»  QR31  =  [zeros(nc,nc) ,zeros(nc,l) ,eye(nc)] ; 

»  save  init  kl  k2  k3  QLll  QRll  QL21  QR21  QL31  qR31 
»  clear  kl  k2  k3  QLll  QRll  QL21  QR21  QL31  QR31  Ac  Be  Cc  Dc  test 

We  now  transform  the  standard  control  problem  into  the  equivalent  decentral¬ 
ized  static  output  feedback  framework  problem  via  the  realization  (4.16). 

>>  A  =  [A,zeros (np,nc) ;zeros(nc,np) ,zeros(nc,nc)] ; 

>>  Bu  =  [zeros(np,nc) ,zeros(np,nc) ,B;eye(nc) ,eye(nc) ,zeros(nc,l)] ; 

»  Bw  =  [Dl;zeros(nc,2)] ; 

»  Cy  «  [zeros(nc,np) .eye(nc) ;C»zeros(l,nc) ;zeros(nc,np) »eye(nc)] ; 

»  Cz  *=  [El  ,zeros(2,nc)]  ; 

>>  Dyu  =  [zeros  (nc  ,nc ), zeros  (nc  ,110) , zeros  (nc ,  1 );  zeros  ( 1  ,nc)  ,  zeros  ( 1  ,nc 

zeros(nc,nc) , zeros (nc .nc) ,zeros(nc, 1)] ; 

>>  Dyw  =  [zeros(nc,2) ;D2;zeros (nc,2)]  ; 

>>  Dzu  =  [zeros (2 .nc) .zeros (2 .nc) ,E2]  ; 

>>  Dzw  ■  EO; 

»  clear  B  C  D  D1  D2  El  E2  EO  nc  np 
>>  who 

Your  variables  are: 

A  Bw  Cz  Dyw  Dzw 

Bu  Cy  Dyu  Dzu 

\\*  ih'W  rntrr  the  remaining  varial*h>  that  ar»>  needful  hv  optgain  (see  Chap- 
!»  r  7»  and  ^ave  them  to  a  file 

>>  k index  ■  [1,1,1]; 

>>  indexkc  •  [4.1.4]; 

>>  indexkr  ■  [4,4,1]; 

>>  N  *  3: 

>>  c t ype  ■  1 ; 

>>  who 

Your  variables  are: 


A 

Cy 

Dyw 

N 

indexkr 

Bu 

Cz 

Dzu 

ctype 

k index 

Bw 

Dyu 

Dzw 

indexkc 

>>  save  datafile 

\W  load  the  file  init. mat  eontaimng  (uir  initial  value  for  the  controller  and 
datafile .  mat  cont  aining  t  In-  (h‘<'rnt  ralized  si  at  ic  output  feedback  problem  data, 
and  >ave  the  combined  data  to  a  singh*  file,  runfile.mat. 
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»  clear 

»  load  datafile 
»  load  init 
»  who 

Your  variables  are: 


A 

Cz 

Dzu 

QL31 

ctype 

k2 

Bu 

Dyu 

N 

QRll 

indexkc 

k3 

Bw 

Dyw 

qLii 

QR21 

indexkr 

k index 

Cy 

Dzu 

QL21 

QR31 

kl 

»  save  runfile 

We  now  clear  the  workspace  and  execute  optgain. 

>>  clear 

>>  [fmin,info,noits]  =  optgain ( 'runfile ’outfile ’ ) 
USING  LINE  SEARCH  ALGORITHM 

ITERATE  FUNCTION  VALUE 


0 

0.4106206586956E+00 

1 

0.4055239785049E+00 

2 

0. 35535104 22910E+00 

3 

0. 35505923 19726E+00 

4 

0.3538665020519E+00 

5 

0 . 3S3246S844085E+00 

ATTEMPTED  STABILITY 

CONSTRAINT  VIOLATION 

6 

0.3523476G21202E*^00 

ATTEMPTED  STABILITY 

CONSTRAINT  VIOLATION 

7 

0.351375379e696E+00 

19 

0.33R5243413508E+00 

20 

0.333t794494535E+00 

21 

0 . 324609 1693388E+00 

22 

0.31947420U990E+00 

23 

0.31689171 00325E+00 

24 

0.3157312171398E+00 

25 

0.3153669369866E+00 

26 

0.3153121011808E+00 

27 

0.3153090709009E+00 
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28  0.3153090381912E+00 

TOTAL  ITERATIONS  =  28 

UNCHND  WARNING  —  INFO  =  1:  PROBABLY  CONVERGED,  GRADIENT  SMALL 
fmin  = 

0.3153 
info  =  ' 

1 

noits  = 

28 

Note  the  warning  ATTEMPTED  STABILITY  CONSTRAINT  VIOLATION  appears  when¬ 
ever  the  optimization  algorithm  was  required  to  reduce  the  step  length  in  order 
to  satisfy  a  constraint.  This  warning  will  appear  occaisonally  during  normal 
optimization.  Similar  warnings  appear  during  scaled  Popov  synthesis. 

We  now  load  the  file  outfile.mat  and  examine  the  optimized  controller 
parameters. 

>>  clear 
>>  load  outfile 
>>  who 

Your  variables  are; 

fmin  info  koptl  kopt2  koptS 

»  koptl 
koptl  “ 

-0.0920 
>>  kopt2 
kopt2  * 

0. 1566 
>>  kopt3 
kopt3  » 

0.  1566 

.\s  a  final  note,  the  function  dsof format  could  have  been  used  to  transform 
the  standard  problem  into  the  e<pn valent  decentralized  static  output  feedback 
framework  problem;  after  entering  the  standard  control  problem  plant  and  gen¬ 
erating  and  saving  an  initial  se't  of  parameters  (kl,  k2,  and  k3),  we  execute 
dsofformat  as  follows 

>>  dsofformat ( Mataf ile ’,lt[l,l] ♦A*B»D1,C»E1,D,D2,E2,E0) ; 

>>  clear 

>>  load  datafile 


>>  who 
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Your  variables  are: 


A 

Cz 

Dztf 

QL31 

indexkc 

Bu 

Dyu 

N 

QRll 

indexkr 

Bv 

Dyw 

QLll 

QR21 

kinder 

Cy 

Dzu 

QL21 

QR31 

Note  that  dsof format  not  only  generates  the  decentralized  static  output  feed¬ 
back  realization,  but  also  the  matrices  QLij  and  QR^j,  and  the  variables  N, 
kinder, indexkc,  and  indexkr. 

4.3  Decentralized  Strictly  Proper  Dynamic  Com- 
pensation 

Lot  ti  and  y  each  be  partitioned  into  two  vector  channels  by  u  -  [  iij  uj 


//  =  [  yj  yj  ]^,  and  rewrite  (4.1)-(4.3)  as 

X  =  Ax  Fi\  u\  D2U2 Diw,  (4.17) 

T/i  =  C\X  F \\tt]  F12U2 (4.18) 

y2  =  Ox  +  »i  +  F22U2  +  7^22^’)  (419) 

z  =  E\X  -f-  £-21*M  4"  ^22^2  4"  Eqw,  (4.20) 

\V*’  fonsidor  a  two-channel,  decentralized,  strictly  proper  dynamic  compensator 

j-,.)  =  .4.-1  Xrl  A  (4.'^1) 

=  rVix.i.  (4.22) 

r^2  =  .4.-'jx,-2  4-  /^<-2i/2«  (4.23) 

U2  =  02^-2  (4-24) 

iHtiirj  K  d<‘tK)te  the  block-diagonal  matrix 

K  =  block-tliag( . 4.  j .  !i,  \ .  t  Vi*  -4c2.  ^’02)-  (4.25) 


It  ran  br  verified  that  nonMiiguiar  Hence,  the  closed-loop  system  con- 

M>ting  of  (4.17)-(4.24)  can  written  d<‘centralized  static  output  feedback 
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with  m  —  V  =  Q,  <f>i  =  <f>2  =  <f>3  —  <1^4  —  —  4^6  —  S’Hd  G(s)  given  by 


A 

0 

0 

0 

0 

0 

0 

52|0|Di  ■ 

0 

0 

0 

1 

; 

0 

0 

0 

0  jo|  0 

0 

0 

0 

0 

0 

0 

1 

I 

0  |0  0 

0 

1 

0 

0 

0 

0 

0 

0 

0  |0|  0 

Cl 

0 

0 

0 

0 

Cii 

0 

0 

^121011921 

0 

/ 

0 

0 

0 

0 

0 

0 

0  |01  0 

0 

0 

/ 

0 

0 

0 

0 

0 

0  jo|  0 

C2 

0 

0 

0 

0 

F21 

0 

0 

F22\^\D22 

0 

0 

/ 

0 

0 

0 

0 

0 

0  |01  0 

0 

0 

0 

_0_ 

0 

0 

0 

0 

0  |oi  0 

Cl 

0 

0 

0 

0 

C21 

0 

0 

B22IOI  Bo 

4.4  Centralized  Strictly  Proper  Dynamic  Com¬ 
pensation  with  Normal  Form  Parameteriza¬ 
tion,  Nominal  Plant 

Consider  a  centralized  strictly  proper  dynamic  compensator  (4.14)  with  the 
(iyn.'iinics  matrix  parameterized  in  normal  form  as 

Ac  —  block*diag(.Aci - -  ^c<5)i 


where 


n, 


i  =  1 . S, 


(4.28) 


ihe  matrices  Be  and  CV  remain  un<'on^^ rained.  Tlie  order  of  the  compen- 
>at()r  i>  tiiiis  26.  Defining 


and  lei t mu  K  denote  the  block-diagonal  matrix 


(4.29) 


K  =  \)\ock-d\hg(<t]  . .  «<5. 6<5,  Be,  C’c),  (4.30) 


it  can  be  verified  that  Lk  i**'  nonsingular.  Hence,  the  closed-loop  system  con¬ 
sisting  of  (4.1)-(4.3)  and  (4.13)-(4.14).  with  the  parameterization  (4.27),  can 
be  written  as  decentralized  static  output  feedback  with  m  =  2<5-|-2,  v  =  S  2, 
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(^1  =  <^2  =  •  •  •  =  =  2,  ^6+1  =  ^*<5+2  =  1,  and  G{s)  given  by 


A 

0 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

5|0pi  ■ 

0 

0 

V^II 

Ki2 

Vll 

-V^12 

••• 

VS2 

Vil 

-v:52 

I 

0  joj  0 

0 

^iT 

0 

0 

0 

0 

...  0 

0 

0 

0. 

0 

0  |0|  0 

0 

0 

0 

0 

0 

...  0 

,  0 

0 

0 

0 

0  |0|  0 

0 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

0  |0|  0 

0 

^iT 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

0  |0|  0 

0 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

0  |01  0 

0 

^(52 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

0  |0|  0 

0 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

0  |0|  0 

0 

V  (5 1 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

0  |0|  0 

c 

0 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

Flo  p2 

0 

/ 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

0  joj  0 

0 

0 

0 

0 

0 

0 

...  0 

0 

0 

Q 

0 

0  |0|  0 

El 

0 

0 

0 

0 

0 

...  0 

0 

0 

0 

0 

F210IF0  . 

(4.31) 


4.5  Centralized  Strictly  Proper  Dynamic  Com¬ 
pensation  with  Uncertain  Plant  Dynamics 

('(insider  the  controller  structure  given  by  (4.13)-(4.14),  and  assume  that 
iiiicert amt\’  in  the  dynamic.s  of  the  plant  is  accounted  for  by  the  model 

X  -  (.1  +  MaA  xNa  i-r  +  Rti  +  Di  w.  (4.32) 


I.etting  K  be  given  by  (4.1*))  and  letting  A  =  Ai  =  A/t,  the  closed-loop  system 
consisting  of  (4.2),  (4.3).  (4.13)  (4  1  l).and  (4.32)  can  be  written  as  decentralized 
static  output  feedback  with  rn  =  i  =  3,  O]  =  Oo  =  (^3  =  1,  p  =  1,  t^’i  =  1.  and 
rr'(s)  given  by 


G(s) 


A 

0 

0 

0 

BIMaIDi  1 

0 

0 

/ 

I 

0  1 

0 

1  0 

0 

I 

0 

0 

0 1 

0 

|0 

r 

0 

0 

0 

El 

0 

IF.. 

0 

I 

0 

0 

0 1 

0 

1  0 

.V4 

0 

0 

0 

0 1 

0 

1  0 

.  £■> 

0 

0 

0 

E21 

0 

IFo  . 

(4.33) 
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4.6  Centralized  Strictly  Proper  Dynamic  Com¬ 
pensation  with  Uncertain  Input  and  Output 
Matrices 

Consider  the  controller  structure  given  by  (4.13)-(4.14),  and  let  the  uncer¬ 
tainty  in  the  input  and  output  matrices  of  the  plant  be  represented  by 

X  —  Ax  (B  ■¥  Mb^bNb)^  +  Diw  (4-34) 


and 

y  =  (C  +  +  Fu  d-  D2W, 

Letting  K  be  defined  by  (4.15)  and  defining  A  as 


A  —  block-diag(AB,  Ac), 


(4.35) 


the  elo.seddoop  system  consisting  of  (4.3).  (4.13),  (4.14),  (4.34),  and  (4.35)  can 
hr  written  as  decentralized  static  output  feedback  with  771  =  ^  =  3,  (l>i  ~  <!>2  — 
Oa  =  1,  p  =  2.  Cl  =  1.  and  6'(s)  given  by 


(4.36) 


4.7  Multiple  Model  Compensation 


r.  .tj^id»  r  a  stri('tly  j)ropf‘r.  <'ciitr;i!i/»'d  »!>  nariih' (  ompensator 

X.  -  . L  -r  -f-  li  y.  (4.37) 

^  Cx  (4.38) 

uhe  h  IN  to  stabih/*  and  t>roMd»'  a.'rrptabh'  performance  for  two  dis- 

tue  1  plarit>.  namely. 

ii  =:  /)nu\  (4.39) 

yj  zr  (\t  f  \  U  D'ilti',  (4.40) 

z\  =  L'l  1 X  +  L'ji  *<  +  Eol  ic,  (4-41) 


and 


.4'jx  -r  B  'jfi  d- 
( *jx  4-  /du  d-  D'n 

El'jX  d-  E'j'jU  d-  E()2«’- 


(4.42) 

(4.43) 

(4.44) 
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Such  multiple  models  may  be  used  to  represent  failure  modes  of  the  nominal 
plant,  or  represent  the  plant  dynamics  at  various  values  of  an  uncertain  param¬ 
eter  [1]. 

Letting  K.  denote  the  partitioned  matrix 

0  0  0  0  0 

0  .4c  0  0  0  0 

0  0  Be  0  0  0 

0  0  0  Be  0  0 

0  0  0  0  Co  0 

0  0  0  0  0  Cc 

the  closed-loop  systems  consisting  of  (4.39)-(4.44)  and  (4.37)-(4.38)  can  be 
written  as  decentralized  static  outj)ut  feedback  with  m  =  6,  v  —  <f>i  —  <p2  — 
03  =  2.  and  Ga(s)  given  by 


(4.45) 


(4.46) 


Chapter  5 

Homotopy  Theory 


5.1  Probability- One  Homotopy  Methods 

Honiotopies  are  a  traditional  part  of  topology,  and  have  found  significant 
application  in  nonlinear  functional  analysis  and  differential  geometry  [11].  Ho- 
iiiotopy  methods  are  globally  convergent,  which  distinguishes  them  from  most 
iterative  methods,  which  are  only  locally  convergent.  The  general  ideaof  homo- 
lopy  methods  is  to  make  a  continuous  transformation  from  an  initial  problem, 
which  can  be  solved  trivially,  to  the  target  f)roblem. 

Following  [12],  the  theoretical  foundation  of  all  probability-one  globally  con- 
vrrgrnt  homotopy  methods  is  givm  in  the  following  different ial  geometry  theo¬ 
rem 

Definition  5.1  Ut  V  C  /C’  and  V  C  IV  (x  stts,  and  Ut  p  :  U  x  [0,  1)  x 
\  ^  fp  /m  a  C*  map  p  ts  said  to  /k  tva\  *rs(d  to  Zi  ro  tf  the  Jacobian  matrix 

I)f’  full  tank  ofi 

Tin‘orein  5.1  If  p(a.\.T)  ts  transifisal  tu  zt  ro.  tinn  for  almost  alia  G  U  the 
voip 

PaiX.Jr)  ~  p{a.  A.x)  (5.1) 

ts  also  transversal  to  zero;  i.e.,  u'tth  prolHihthty  one  the  Jacobian  matrix  Dpa(X,  x) 
has  full  rank  on 

The  recipe  for  constructing  a  globally  convergent  homotopy  algorithm  to 
solve  the  nonlinear  system  of  e(|uations 

fU)  =  0.  (5.2) 

wher('  f  \  is  a  f’*  map.  is  as  follows:  For  an  open  set  U  C  R"^ 

construct  a  C“  homotopy  ma[)  p  :  T  x  [0,  1)  x  7?^  — R^  such  that 
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1)  p{a,X,x)  is  transversal  to  zero, 

2)  pa(0,  a:)  =  p(a,  0,  z)  =  0  is  trivial  to  solve  and  has  a  unique  solution  xq, 

3)  pa(l,a;)  =  f(x), 

4)  Pa  is  bounded. 

Then  for  almost  all  a  G  there  exists  a  zero  curve  7  of  pa,  along  which  the 
Jacobian  matrix  J9pa  has  rank  p,  emanating  from  (0,xo)  and  reaching  a  zero  x 
of  /  at  A  =  1.  This  zero  curve  7  does  not  intersect  itself,  is  disjoint  from  any 
other  zeros  of  pa,  and  has  finite  arc  length  in  every  compact  subset  of  [0, 1)  x 
Furthermore,  if  Df{x)  is  nonsingular,  then  7  has  finite  arc  length.  The  general 
idea  of  the  algorithm  is  to  follow  the  zero  curve  7  emanating  from  (0,  xq)  until 
a  zero  x  of  f(x)  is  reached  (at  A  =1). 

The  zero  curve  7  is  tracked  by  the  normal  flow  algorithm  [12],  a  predictor- 
corrector  sceme.  In  the  predictor  phase,  the  next  point  is  produced  using  Her- 
inite  cubic  interpolation.  Starting  at  the  predicted  point,  the  corrector  iteration 
involve.s  computing  (implicitly)  the  Moore-Penrose  pseudo-inverse  of  the  Jaco- 
luan  matrix  at  each  point.  The  most  complex  part  of  the  homotopy  algorithm 
is  the  computation  of  the  tangent  vectors  to  which  involves  the  computation 
of  tlie  kernel  of  the  p  x  (p  +  1)  Jacobian  matrix  Dpa-  The  kernel  is  found  by 
computing  a  QR  factorization  of  Dp^.  and  then  using  back  substitution.  This 
strategy  is  implemented  in  the  mathematical  software  package  HOMPACK  [14], 
which  was  used  for  the  curve  tracking  here. 

1  wo  different  homotopy  maps  are  used  for  solving  the  optimal  projection 
f‘(|uations.  When  the  initial  problem,  5(x;a)  =  0,  can  be  solved,  then  the 
luunotopy  map  is  [13] 

pa(A,  x)  =  F(u.  A.  x)  =  A/(j-)  4*  ( 1  “  A)<7(x;  a),  (■^•^) 

wIm  fix)  =  0  is  the  final  proi>iem.  and  a  is  a  parameter  vector  used  in  defining 
t  Ic  function 

Who’ll  the  initial  problem  is  not  solvanl  exactly,  i.e..  //(xq;^)  ^  0,  then  the 
m  ip  a  Newton  homotopy  [9] 

pa  (A,  x)  =  F{h.  A.  x)  -  ( 1  -  A)/‘  {b.  0,  xo),  (0-4) 

uh-P-  (I  -  (/xxo).  For  A  =:  0.  /*a(l).xo)  -  F(6,U,xo)  =  0,  and  for  A  =  1, 
/.pl.x.»^F(6,l,x)  =  /(x)^() 

When  the  map  (5.4)  is  ised,  the  <'<piations  considered  for  0  <  A  <  1  are  not 
the  optimal  projection  equations  whenever  F(6,0,xo)  =  p(xo;t)  ^  0.  Hence,  a 
goal  in  constructing  the  initial  system  and  the  starting  point  may  be  to  minimize 


Chapter  6 

Formulation  (Homotopy) 


6.1  The  order  reduction  problem 

The  'Ho-optimal  model  order  reduction  problem  is  that  of  approximating  a 
higher  order  dynamical  system  by  one  of  lower  order  so  that  a  quadratic  model 
rediK'tion  criterion  is  minimized. 

Idle  problem  can  be  formulated  a.s:  given  the  asymptotically  stable,  control¬ 
lable.  observable,  time  invariant,  continuous  time  system 

i(t)  Ar{t)  +  Bu(t),  (6.1) 

y{i)  =  rr(l).  (6.2) 

wlicr*'  A  e  B  e  TV'"'”' R‘”‘.  tlio  goal  is  to  find  a  reduced  order 

iii-.d.-l 

i„.(0  =  A,.,T„An  +  B.„u(t).  (6.3) 

(6.4) 

wl,<n  B,„  e  R"-’'”’.  (',n  €  R'’"'”'.  <  n,  which  minimizes 

tip  funrtion 

J ( .-l,,. .  B,„ .  C,n )  H  [(.'/-  i/.M  V  R(y-  Vm )] ,  (6.5) 

wher<‘  the  input  u(<)  is  white  nois<‘  with  symmetric  and  positive  definite  intensity 
r  and  /?  is  a  symmetric  and  positive  definite  weighting  matrix. 

6.2  The  combined  H  model  reduction  prob¬ 
lem 

In  practice,  to  simplify  a  plant  for  control  design  or  to  simplify  a  controller 
for  ease  of  implementation,  a  role  must  be  taken  into  account,  i.e.,  the  order 
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reduction  approach  should  approximate  the  system  frequency  response  to  the 
greatest  extent  possible. 

The  problem  is  formulated  as:  given  the  controllable  and  observable,  time 
invariant,  continuous  time  system 


x{t)  =  Ax(t)  +  BDu{t),  ■  (6.6) 

y(<)  =  Cx(t),  (6.7) 

where  t  £  [0,oo),  A  £  is  asymptotically  stable,  B  £  C  £  R‘", 

D  £  R'"’‘f’  (m  <  p)  and  the  input  Du(t)  is  white  noise  with  symmetric  and 
positive  definite  intensity  V  =  DD'^  ,  find  a  rim-th  order  model  (n^  <  «) 

im{t)  =  Am  +  Bm  Du{t),  (6.8) 

ym(<)  =  CmXmit),  (6-9) 


where  €  R"™ Bm  £  R"’"’‘'".  C,,,  £  R-'"”' .  which  satisfies  the  following 
criteria: 

(i)  Arn  is  asymptotically  stable; 

(ii)  the  transfer  function  of  the  reduced  order  model  lies  within  7  of  the  transfer 
function  of  the  full  order  model  in  the  11^  norm,  i.e., 

\\H[s)  -  //„,(s)|U  <  7  (6.10) 

where 

//(.*»)  =  EC  {sIri  —  -4)  ^  BD.  =  EC rn{s Im  BmD^  (6.11) 

-  >  0  IS  a  given  constant.  E  £  Ci  >  /)  is  a  given  constant  matrix;  and 

(m)  thr  //‘  model  reduction  criterion 

J{Am-  Brn.Cm)  =  hill  [(i/  -  y,„)^/f(.v  -  y,,,)]  (6.12) 

i>  fniniiiHzetl.  where  S  is  the  expected  value  and  H  —  E^ E  is  a  symmetric  and 
[K.'Nitivf  definite  weighting  matrix. 

6.3  The  LQG  controller  synthesis  with  an  Tioo 
performance  bound 

Tlic  'Hi/'Hca  mixed-norm  controller  synthesis  problem  provides  the  means 
for  simultaneously  addressing  Hi  and  H^  performance  objectives.  In  practice 
such  controllers  provide  both  nominal  performance  (viaTfo )  and  robust  stability 
(via  H^  )■  Hence  mixed-norm  synthesis  provides  a  technique  for  trading  off 
performance  and  robustness,  a  fundamental  objective  in  control  design.  (It 
should  be  noted  that  Hn  controller  synthesis  is  a  special  case  of  mixed-norm 
controller  synthesis,  with  the  H^  bound  set  to  oo). 
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The  LQG  controller  synthesis  problem  with  an  H°°  performance  bound  can 
be  stated  eis:  given  the  n-th  order  stabilizable  and  detectable  plant 

x{t)  =  Ax(t)  +  Bu{t)  +  Diw(t),  (6.13) 

y(t)  =  Cxit)  +  D2w{t),  (6.14) 

where  A  G  B  €  C  G  R'’'",  Di  G  R’”'p,  D2  G  R'^'p,  DiD^  =  0, 

and  w{t)  is  p-dimensional  white  noise,  find  a  Uc-th  order  dynamic  compensator 

i,{<)  =  Acx,(t}  +  Bcy{t),  (6.15) 

u{i)  =  CcXc(<).  (6-16) 

where  .4c  G  Be  G  R'’^'^'.  Ce  G  R"”""',  and  Uc  <  n,  which  satisfies  the 

following  criteria: 

(i)  the  closed-loop  system  (6.13)  -  (6.16)  is  asymptotically  stable,  i.e., 

= ( ,ic  T ) 


i.-s  a.s\  in[)tot ically  stable: 

(ii)  tli<‘  7-*^  X  p  closed-loop  transfer  fnnrlion  from  U‘(0  to  E:c^x{t)  £'2cow(0t 


//(.) 

where 

/:  .  =  {  E.^Cr  )  (  /:ix 


vinsf\  tie  rttiistraint 


( ••«/,■,  -  .4)-' 

D. 

(6.18) 

€  R' 

GR*'-"'".  Ele^E2c^  =  0), 

(6.19) 

h  = 

11  It, 

(6.20) 

/>=( 

Ih  ^ 

B  e.  / 

(6.21) 

i7/i> 

i  ,  - 

(6.22) 

ule  r‘  “  '  U  1''  a  e,ivf  n  ('(»n>!ant  and 

niM  lie  [>«Th>rnian<'e  funet i< inal 


J(Ac^BcX\) 


^liin  (  /  )/{i  j(/)  + 


(6.23) 


i>  miiMinized,  w’here  €  is  the  exp^^Med  value,  R\  =  EJ E\  G  and  R2  — 

El  E’>  €  (E,  G  E'j  €  Ej  E^  =  0  )  are  respectively  sym- 

nietrie  positive  seniidefinite  arul  symttietrie  positive  definite  weighting  matrices. 
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6.4  Functions 

The  thirteen  functions  listed  here  are  classified  by  their  purpose.  (They  are 
described  in  detail,  with  examples  of  their  use,  in  Chapter  7). 

morh2inf,  inorh21y,  morli2over,  and  morh2op  are  for  the  ?f2-optimal 
model  order  reduction  problem.  morh2hiinf,  morh2hily,  and  morh2hiover 
are  for  the  combined  ?f2/7fco  model  order  reduction  problem.  For  the  reduced 
order  LQG  ('H2-optimal  )  problem  use  rlqgly,  rlqginf,  or  rlqgover.  To  solve 
the  full-order  LQG  ('H2-optimal )  controller  synthesis  problem  with  an  'Hoc  norm 
bound,  use  flqgly,  flqginf,  or  flqgover. 

6.5  Final  Note 


[fi]. 


For  a  complete  explanation  of  the  homotopy  algorithms  described  above,  see 


Chapter  7 

Program  Descriptions 
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dfolqg,  drolqg 


Purpose 

To  synthesize  full-  and  reduced-order  discrete-time  ?{2-optinial  dy¬ 
namic  compensators. 

Synopsis 

[Ac,Bc,Cc,Dc,cost3  =  dfolqg(A,Bu,Bw,Cy ,Cz,Dyu,Dyw,Dzu,Dzw) 
CAc.Bc.Cc.Dc.test]  =  drolqg(A.Bu,Bw,Cy,Cz,Dyu,DyM,Dzu,Dzw,nc) 

Description 

dfolqg  synthesizes  discrete-time  T^o'Optimal  compensators  based  on 
routines  in  the  Matlab  Control  Systems  Toolbox.  Given  an 
nth-order  two- vector  imput,  two  vector  output  plant 

x(A* -h  1)  =  Ax(k)  Duu{^) 

y{k)  =  Cyx{k)  Dyuu(k)  -h  Dyuiw{k), 
z{k)  =  C:X(k)  +  ll(k)  -h  D,u>w(k), 

dfolqg  designs  a  full-order  'W2-<^Ptimal  strictly-proper  dynamic  com¬ 
pensator,  of  the  form 

Xc{k  -f  1)  -  .1  )  -h  Bcy(k) 

u{k)  =  i'rxAk)-^  IXy(k). 

S<>\r  that  the  realization  A,  .Br  (\  Bc  returned  by  dfolqg  is  not 
tie*  same  as  that  produ('ed  b\  dreg  m  the  CONTROL  SYSTEMS 
loot. BOX.  which  usf»s  llie  rurn  ut  estimator  of  [5],  while  dfolqg 
the  ;urd/r/or estimator  the  (‘ONTROL  SYSTEMS  lOOLBOX 

t  vr  ^  (iuide  for  detail  [7]  cost  i.s  the  value  of  the  discrete-time 
H\:  norm  for  the  closeddoop  transfer  function 

drolqg  uses  a  similarity  tran>f()rmation  to  balance  the  realization 
of  the  controller,  and  then  eliminates  n  —  of  the  least  observ¬ 
able  and  controllable  state's  in  this  realization  in  order  to  obtain 
a  reduced-order  approximation  of  the  full-order  compensator.  Bal¬ 
anced  realizations  do  not  exist  for  unstable  controllers,  and  drolqg 
will  fail  with  an  error  message  if  dfolqg  returns  an  unstable  compen- 
.sator.  test  is  a  flag  which  signals  whether  the  balanced  truncation 
returned  by  drolqg  stabilizes  the  closed-loop  system. 


test  =  0  — f  closed-loop  stable, 

test  =  1  — ^  closed-loop  unstable. 
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Notes 

See  dlqr  and  dlqe  from  CONTROL  SYSTEMS  TOOLBOX. 

Examples 

»  A  =  [1.0000  0.0001  0  0;  -0.000025  1.0000  0  0;  0  0  1.0  0.0001; 
0  0  -0.0004  1.0000] ; 

»  B  =  [0,  0;  0.00005,  0.00005;  0,  0;  0.0001,  0.00020]; 

»  D1  =  [B,  zeros(4,2)]; 

»  C  =  [0  1  0  -0.5;0  1  0  -1.0] ; 

»  El  =  [C;  0  0  0  0;  0  0  0  0]; 

»  D  =  zeros(2,2) ; 

»  D2  =  [0,  0,  0.1,  0;  0,  0,  0,  0.1]; 

»  E2  =  [0 ,  0 ;  0 ,  0 ;  0.1,  0 ;  0 ,  0.1]; 

>>  EO  =  zeros(4,4); 

>>  who 

Your  variables  are: 

A  C  D1  EO  E2 

B  D  D2  El 

»  tAc»Bc,Cc,Dc,cost]  «  dfolqg(A,B*Dl ,C,E1 ,D,D2,E2,E0) 

Ac  « 


1.0000 

-0.0003 

0 

0.0002 

0.0001 

1.0004 

-0.0004 

-0.0009 

0 

-0.0004 

1.0000 

0.0003 

0.0002 

0.0037 

-0.0012 

0.9958 

0.0003 

0.0001 

0.0004 

-0.0006 

0.0004 

-0.0000 

-0.0001 

-0.0021 

m 

1.5604 

-7.6711 

-8.2375 

-0.4167 

0.4664 

11.4192 

0.2584 

-10.3769 

Dc  « 

0  0 

0  0 

cost  = 

2.3624e-05 
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>>  [Ac,Bc,Cc,Dc,test]  =  drolqg (A,B ,D1 ,C,E1 ,D,D2 ,E2,E0 ,2) 
Closed-Loop  System  UNSTABLE 
Ac  = 


1.0000 

0.0000 

-0.0000 

1.0000 

0.0005 

0.0068 

-0.0341 

0.0177 

0.0005 

0.0341 

0.0068 

-0.0177 

Dc  - 

0  0 

0  0 

test  = 

1 
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dsofformat 


Purpose 

To  transform  standard  control  problem  into  the  equivalent  decen¬ 
tralized  static  output  feedback  framework  problem. 

Synopsis 

dsofformat (filename ,ncvec,noivec,A,Bu,Bw,Cy , 

Cz , Dyu , Dy w , Dzu , Dzw) 

dsof format (filename ,ncvec ,noivec ,  A,Bu,Bw,Cy , 

Cz , Dyu , Dy w , Dzu , Dzw , Ma , Na . Mb , Nb , Me , Nc ) 

Description 

dsofformat  reformulates  a  A*-(dianiiel  fixed-structure  control  prob¬ 
lems  into  the  equivalent  de<-ent ralized  static  output  feedback  frame¬ 
work  problem. 

nevee  is  a  A*  X  1  vector  whose  i-\  h  element  is  the  order  of  the  dynamics 
for  the  ;-th  processor.  Stali('  gain  compensators  are  indicated  by  a 
dimension  of  0. 

noivec  is  a  kx2  matrix.  The  ?*th  elf'ment  of  the  first  column  con¬ 
tains  the  dimension  of  mejisureuient  signal  for  the  /-th  channel,  and 
tie"  /-th  element  of  the  M*c(,)nd  <'o!nmn  contains  the  dimension  of  the 
aetualcir  signal  for  tie*  i-lli  <  lianre  1 

If  t  h<*  arguments  Ma.  Na.  ♦  tc  .  ar»*  iie  ludt  d.  dsofformat  will  generate 
th'‘  matrices  required  to  repre>#  rit  th#*  fr>l!<aving  plant  uncertainty  in 
th'  d»<**ntralizf*<l  state*  output  frvdbaek  framework 


;  H,,  -r 

(  tj  -  M.  Ao  .V, 

;  iKu  J 

I  le*  output  of  dsofformat  conM>ts  of  the  matrices  of  the  decentral- 
iz^*<l  static  output  feedback  reahzaticui  (2.2),  the  matrices  Qhij 
Qiuj  which  define  the  structure  of  the  controller  (2.14),  and  the  rep¬ 
etition  and  sizing  variables  N  kindex.  indexkc,  and  indexkr  (see 
optgain  for  details  on  thesr  variabh*s).  1  hese  output  variables  are 
saved  to  the  file  fihtiamf  .mat 

Notes 

set*  also  optgain 
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Examples 

First  we  will  demonstrate  setting  up  single-channel,  reduced-order  central¬ 
ized  strictly-proper  dynamic  compensator  problem: 

»  clear 

»  A  =  [zeros (5 , 1) ,eye (5) ; -1  -2  -24  -12  -24  -4]; 

»  B  =  [zeros  (5 , 1) ; 1] ; 

»  C  =  [1  0  0  0  0  0]; 

»  D  =  0; 

»  D1  =  [B,zeros(6, 1)] ; 

»  D2  =  [0,1]  ; 

»  El  =  [C;zeros(l,6)] ; 

»  E2  =  [0;1]; 

>>  EO  =  zeros (2,2); 

>>  who 

Your  variables  are: 

A  C  D1  EO  E2 

B  D  D2  El 

>>  dsof format ( *data \4,[l,l] ,A,B,D1,C,E1,D,D2,E2,E0) ; 

>>  clear 
>>  load  data 
>>  who 

Your  variables  are: 


A 

Cz 

Dzw 

QL31 

indexkc 

Bu 

Dyu 

N 

QRll 

indexkr 

Bw 

Dyw 

QLll 

QR21 

k index 

cy 

>> 

Dzu 

QL21 

QR31 

I  Im  containdf'd  iii  data  . mat  an-  tlif'  transformed  matrices  correspond- 

MIC  t(.  tlic  n-aliziition  (2.‘J)  I  li-  iii.tiri.  o  QLij  and  QRy  correspond  to  (2.14). 

A  s*:'cond  example  shows  an  exampi*'  of  .s<*tting  up  a  two-channel  decentral- 
iz<'d  dynamic  compensator  syntlu'sis  problem  in  the  decentralized  static  output 
fi'ctlhack  framework,  including  real-parameter  uncertainty  in  the  input  matrix. 

>>  clear 

»  A  =  [0,  l:-3  -4] : 

»  B  »  [0,  0:-l,  -.31: 

»  Dl  =  [35.  0,  0;  -61.  0.  0] ; 

»  C  =  [2.  1:3.  1]  : 
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»  El  =  [52.1950,  8.9440;  0,  0;  0,  0]; 

»  D  =  zeros (2, 2) ; 

»  D2  =  [0,  1,  0;  0,  0,  1]  ; 

»  E2  =  [0,  0;  1,  0;  0,  1]  ; 

»  EO  =  zeros(3,3) ; 

»  Mb  =  B; 

»  Nb  =  eye  (2) ; 

>>  dsof  format  ( ’datafile  ’ ,  [2  ;2]  »  [1 » 1 ;  1 , 1]  ,  A  ,B,D1  ,C,E1  ,D ,  D2  ,E2  ,E0 ,  []  ,  []  ,  Mb,  Nb,  []  ,  []  ) ; 
>>  clear 

»  load  datafile 
>>  who 

Your  variables  are: 

A  Cz  Dyw  QL31  QR31  k index 

Bd  Ded  Dzu  QL41  QR41 

Bu  Deu  Dzw  QL51  QR51 

Bw  Dew  N  QL61  QR61 

Ce  Dyd  QLll  QRll  indexkc 

Cy  Dyu  QL21  QR21  indexkr 

» 

For  more  examples  on  using  dsofformat,  see  Chapter  4. 
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fiqgly,  flqginf,  fiqgover 


Purpose 

Find  the  full-order  LQG  compensator  with  an  Tioo  norm  bound. 
Synopsis 

[Ac,  Be,  Cc,  cost]  =  flqgly(A,  B,  C,  D,  gammaO,  gamma,  El,  E2,  Eli, 
E2i,  Dl,  D2) 

[Ac,  Be,  Cc,  cost]  =  flqginf(A,  B,  C.  D,  gammaO,  gamma,  El,  E2, 

Eli.  E2i,  Dl,  D2  ) 

[Ac,  Be,  Cc,  cost]  =  fiqgover  (A,  B,  C,  D,  gammaO,  gamma,  El,  E2, 
Eli,  E2i,  Dl,  D  2) 

Description 

For  a  given  order  linear  plant  with  open-loop  state  space  real¬ 
ization  given  by 


i(0  = 

.4jr(0  +  Du{1)  +  Diu'it) 

(7.2) 

y(t)  = 

Cx(l)  +  DMt) 

(7.3) 

:{t)  = 

£,/•(/)  +  E2«(<) 

(7.4) 

ec(0  = 

(7.5) 

fiqgly,  flqginf.  and  fiqgover  calculate'  Ac,  Be,  and  state 

vpat'e*  realization  for  the  fullorde'r  LQ(i  ( Wj  optimal)  compensator 
winch  yields  a  closed-loop  syst#  in  with  H  ^  norm  bounded  by  gamma, 
l  ie*  closed-loop  Tij  cost  is  given  by  cost  gammaO  (-/o)  is  the  initial 
-  and  should  always  be  greater  than  gcunma. 

lie*  rf'sulting  compensator  from  tt  fhiginf  is  in  the  input  normal 
Hiccati  form  while  tliat  from  fiqgly  is  m  by  s  form. 

F  X  a  m  p  les 

>>  A*2eros(8); 

>>  B»zero8(8, 1) ; 

»  C»zeros(l ,8) ; 

»  D=0; 

»  A(l:8,l)  =  [-0.161:  -6.004;  -0.5822;  -9.9835;  -0.4073;  -3.982; 

0; 

0;]; 

>>  for  i  ==1:7  A(i,i+1)  *  1;  end 
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»  B(l:8,l)  =  [0;  0;  0.0064;  0.00235;  0.0713;  1 .0002;0. 1045;  0.9955;]; 

»  C(l,l)  =  1-0; 

»  El  =  zeros (2, 8) ; 

»  Eld, 1:8)  =  0.001  *  [0  0  0  0  0.55  11  1.32  18]; 

»  Eli  =  El; 

»  E2  =  [0 ;  1]  ; 

»  E2i  =  [0;0]  ; 

»  D1  =  zeros(8 ,2) ; 

»  Dl(l:8,l)  =  B; 

»  D2  =[0  1]; 


^  [Ac, Be, Cc, cost]  = 

flqgly(A,B, 

,C,D,1.0€3, 

,2.0,E1,E2, 

,Eli,E2i,Dr 

,D2) 

Columns  1 

through  7 

0.0000 

1.0000 

-0.0000 

0.0000 

-0.0000 

-0.0000 

0.0000 

-3.4454 

-0.1189 

-0.0000 

0.0000 

-0.0000 

0.0000 

0.0000 

-0.0000 

-0.0000 

0.0000 

1.0000 

0.0000 

-0.0000 

-0.0000 

0.0000 

0.0000 

-1.9890 

-0. 1688 

0.0000 

0.0000 

-0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-0.0000 

1.0000 

0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.6728 

-0.1541 

-0.0000 

0.0000 

0.0000 

0 

-0.0000 

0.0000 

-0.0000 

0.0000 

-0.0000 

0.0000 

0.0000 

0.0000 

-0.0000 

0.0000 

-0.3032 

Column  8 
0.0000 
0.0000 
-0.0000 
-0.0000 
0.0000 
-0.0000 
1.0000 
-0.8793 
Be  - 

0.0012 

0.0048 

0.0061 

0.0094 

0.0068 

0.0143 
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-0.1278 

0.1000 

Cc  = 

Columns  1  through  7 

1.0000  0.0000  1.0000  0.0000  1.0000  -0.0000 
Column  8 
0.0000 
cost  = 

0.1434 

Algorithm 

The  algorithms  for  flqgly,  flqginf ,  and  flqgover  are  described 
in  Chapters  14,  15,  16,  17,  and  18  of  [6]. 

See  Also 


1.0000 


rlqgly.  rlqginf ,  rlqgover 
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folqg,  rolqg 


Purpose 

To  synthesize  full  and  reduced-order  continuous-time  ?{2-optimal  dy¬ 
namic  compensators. 

Synopsis 

[Ac,Bc,Cc,Dc,cost] 

[Ac,Bc,Cc,Dc,test] 

Description 

folqg  synthesizes  continuous-lime  H2-optimal  compensators  based 
on  routines  in  the  Matlab  Control  Systems  Toolbox,  Given 
an  rah-order  two-vector  impul,  two- vector  output  plant 

i(t)  =  Ax(t)-^[Ku(t)-^Buw(t). 

y(t)  =  CyX(t)  -h  Dv„U(/)  +  DyU'W(f), 

z(t)  =  + 

folqg  designs  a  full-order  strictly-proper  dynamic  com- 

p»‘nsator.  of  the  form 

r,(t)  =  A,jrAt]^fLy(t) 

u(t)  -  (\  x  {t\  ^  !),y(f). 

In  >tatulHr(l  LQG  tli(‘ory.  D  ~  0  cost  i>  the  value  of  tlie  'H-j  norm 
\‘  T  tie*  <'l(»>ed*loo|)  tran^hT  fun*  t ion 

rolqg  u>»*>  a  .similaritx  t ran>forfnai n ui  to  balance  tin'  realization 
. -f  tin  controlhr.  and  tln-n  ♦jnmnato  n  —  u,  of  the  h'ast  observ- 
abl'  and  controiiable  statt^^  in  thi.s  realization  in  order  to  obtain 
a  r**du*  ed'order  approxnnat n »n  of  tin-  full-order  comi)ensator.  Bal- 
ain  ed  realizations  do  in^t  exi>t  for  unstable  controllers,  and  rolqg 
will  fail  with  an  error  iiu^sage  if  folqg  returns  an  unstable  compen¬ 
sator.  test  is  a  flag  which  signals  whether  the  balanced  truncation 
returned  by  rolqg  stabilizes  the  clos#^i-loop  system. 

test  =  0  closetl-loop  stable, 

test  I  — >  closed-loop  unstable. 


=  folqg(A,Bu,Bw,Cy,Cz,Dyu,Dyw,Dzu,Dzw) 

=  rolqg(A,Bu,Bw,Cy,Cz,Dyu,Dyw,Dzu,Dzw,nc) 


Notes 

Sef‘  iqr,  Iqe,  and  reg  from  C’ontrol  Systems  Toolbox. 
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Examples 


» 

A  =  [zeros(5, 1) ,eye(5) ;-l 

-2  -24  -12  -24  -4] ; 

» 

B  =  [zeros (5, 1) ;  1]  ; 

» 

C  =  [1  0  0  0  0  0]  : 

» 

D  =  0; 

» 

D1  =  [B, zeros (6,1)] 

t 

» 

D2  =  [0,1]; 

» 

El  =  [C;zeros(l,6)] 

; 

» 

E2  =  [0;1]; 

» 

EO  =  zeros(2,2) ; 

» 

[Ac,Bc,Cc,Dc,cost] 

=  folqg(A.B,Dl,C,El,D,D2,E2,E0) 

Ac 

= 

-0.1480  1.0000 

0 

0  0  0 

-0.0110  0 

1.0000 

0  0  0 

0.0103  0 

0 

1 . 0000  0  0 

0.0016  0 

0 

0  1.0000  0 

-0.0041  0 

0 

0  0  1.0000 

-1.4149  -5.4340 

-25.9960 

-15.5854  -24.6029  -4.1480 

Be 

s 

0.1480 

0.0110 

-0.0103 

-0.0016 

0.0041 

0.0007 

Cc 

m 

-0.4142  -3.4340 

-1.9960 

-3.5854  -0.6029  -0.1480 

Dc 

m 

0 

cost  ■ 

0.2822 

» 

[Ac.Bc.Cc.Dc.test] 

■  rolqgiA 

.B.D1.C,E1,D,D2,E2,E0,4) 

Ac 

» 

-0.0157  -0.3137 

0.0046 

-0.0016 

0.3137  -0.3487 

0.027S 

-0.0093 

-0.0046  0.0275 

-0.2634 

1.0135 

-0.0016  0.0093 

-1.0135 

-0.0360 
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Be  = 

-0.0763 
0.2845 
-0.0112 
-0.0038 
Cc  = 

-0.0763  -0.2845  0.0112 

Dc  = 

0 

test  - 
0 


-0.0038 
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morh2inf,  morh21y,  morh2over,  morh2op 


Purpose 

Find  the  optimal  reduced-order  model  of  a  linear  system. 
Synopsis 

[Am,  Bm,  Cm,  cost]  =  morh2inf (A,  B,  C,  nm) 

[Am,  Bm,  Cm,  cost]  =  morh21y(A,  B,  C,  nm) 

[Am,  Bm,  Cm,  cost]  =  morh2over(A,  B,  C,  nm) 

[Am,  Bm,  Cm,  cost]  =  morh2op(A,  B,  C,  nm,  meth,  init,  cl,  c2) 

Description 

For  a  given  linear  system  with  state  space  representaion  .4,  B,  and  C, 
morh2inf ,  morh21y,  mcrh2over,  and  morh2op  return  the  Tin  optimal 
r^xluced-order  model  .4,,,,  /im-  and  Cm  of  dimension  nm  with  Uo  cost 
cost.  The  result  from  n.orh2inf  is  in  the  input  normal  form  while 
that  from  morh21y  is  in  Ly's  form, 

in  morh2op,  meth  and  init  denote  the  strategy  and  the  method  of 
initialization  respectively  [15].  and  cl  and  c2  define  the  initial  A  by 
4o  =  -cl  /  +  c2.4. 

F  X  am  pies 

>>  A  •  C-IO  1  0;  -5  0  1;  -1  0  0]; 

>>  B  ■  [0  ;  1  ;  1]  ; 

>>  C  *  [1  0  0]  ; 

>  >  nm  ■  2  ; 

>>  [Aa.  Bm.  Cm,  cost  ]  *  aorh2inf(A.  B,  C,  nm) 

Am  « 

“0.1 39T  -0.1006 

0.6010  -0.4482 

Ba  • 

-0.5285 

0.9468 

Cm  * 

-0.3204  -0.0961 

cost  = 

3.2902e-04 

>>  [Am,  Bm,  Cm,  cost  3  *  aorh21y(A,  B,  C,  nm) 
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Am  = 

0  1.0000 

-0.1231  -0.5878 

Bm  = 

0.0784 

0.0782 

Cm  =  . 

1 . 0000  0 . 0000 

cost  = 

3.2902e-04 

»  A  =  [-1  3  0;  -1  -1  1;  4  -5  -4] ; 

>>  B  =[  -2;  2;  4] ; 

»  C  =  [10  0]; 

>>  nm  =  1 ; 

>>  [Am,  Bm,  Cm,  cost  ]  =  morh2op(A,  B,  C,  nm,  2,  2,  10.0,  0.0) 

Am  = 

-10.4365 
Bm  * 

-1.5972 
Cm  = 

1.5972 
cost  * 

1 .6882 

A  Iftorit  h  III 

Jhc  algorithms  for  morh2inf ,  morh21y.  and  morh2over  are  described 
in  C'hapters  2.  4,  and  0  respect  i\>‘ly  of  [(>].  The  algorithm  for  morh2op 
ran  be  found  in  [la]. 

AUo 

morh2hiinf,  morh2hily.  morh2hiover 
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morh2hiinf,  morh2hily,  morh2hiover 


Purpose 

Find  the  combined  U^l^oo  reduced-order  model  of  a  linear  system. 

Sy  nopsis 

[Am,  Bm,  Cm,  cost]  =  morh2hiinf (A,  B,  C,  nm,  gamma) 

[Am,  Bm,  Cm,  cost]  =  morh2hily(A,  B,  C,  nm,  gamma) 

[Am,  Bm,  Cm,  cost]  =  morh2hiover (A,  B,  C,  nm,  gamma) 

Description 

For  a  given  linear  plant  with  state  space  representaion  A,  B,  and 
C\  morh2hiinf ,  morh2hily,  and  morh2hiover  return  the  combined 
'Hi /'Hoc  reduced-order  model  Bm.  and  Cm  of  dimension  nm 
witli  H2  cost  cost.  The  triple  (Am .  Bm.  Cm)  returned  from inorh2hi inf 
is  in  the  input  normal  form  while  that  from  morh2hily  is  in  Ly’s 
form. 

Examples 

>>  A  **  zeros(lO) ; 

>>  B  =  zeros ( 10 ,  1) ; 

>>  C  *  zerosd  ,  10) ; 

»  A(l,l:10)  =  [-10  -45  -120  -210  -252  -210  -120  -45  -10  -1]; 

>>  for  1*1:9  A(i+l,i)*  1.0;  end 
>>  B(l,l)  *  1.0; 

>>  C(1,10)  =  1.0; 

>>  nm  ■  4; 

>>  gamma  *  1.0; 

>>  [Am,  Bm,  Cm,  cost  ]*  morh2hi inf  (A ,  B,  C,  nm,  gamma) 

Am  ■ 


-0.0273 

-0.1286 

-0.0274 

0.0124 

0.2376 

-0.1087 

-0.1936 

0.0397 

-0.1352 

0.5178 

-0.2416 

0.2322 

-0.2412 

0.4166 

-0.9124 

-0.4787 

Bm  = 

0.2338 

-0.4663 
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0.6951 

0.9785 

Cm  = 

0.1897  0.2047  0.1142  -0.0409 

cost  = 

1.2868e-04 

»  [Am,  Bm,  Cm,  cost  ]=  morh2hiover (A ,  B,  C,  nm,  gamma) 
Am  = 


-0.0308 

-0.1739 

-0.0642 

0.0566 

0.1739 

-0.1200 

-0.3132 

0.1338 

-0.0642 

0.3132 

-0.2566 

0.4553 

-0.0566 

0.1338 

-0.4553 

-0.4489 

Bm  « 

0.2148 
-0.3180 
0.2916 
0 . 2044 

Cm  = 

0.2148  0.3180  0.2916  -0.2044 

cost  ■ 

1.2868e-04 

A  igorit  h  tn 

Wi*  algorithms  for  morh2hiinf ,  morh2hily,  and  inorh2hiover  are 
d»*>' ril»^‘d  in  C  hapters  Ih  and  10  r^^spert ively  of  [()]. 

s.,  Als«. 

inorh2op,  .morh2inf,  morh21y.  morh2over 
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optgain 


Purpose 

To  optimize  fixed-structure  controllers  using  a  quasi-Newton  gradi¬ 
ent  search  algorithm. 

Synopsis 

[fmin,  info,  noit]  =  optgain(inf ile,  outfile) 

[fmin,  info,  noit]  =  optgain(inf ile ,  outfile,  outon,  diag, 
maxit,  method,  tolfac) 

Description 

optgain  optimizes  fixed-structure  controllers  using  a  quasi-newton 
gradient  search  algorithm.  The  data  defining  the  problem  to  be 
solved  is  contained  in  the  file  ttijilt  and  must  include  the  fol¬ 

lowing: 

•  A,Bu,Bw,Cy ,Cz,Dyu,Dyw,Dzu,Dzw,  the  matrices  of  the  decen¬ 
tralized  static  output  feedback  framework  realization  (2.2). 

•  QLij  ,QRij,  i  =  1 . ^  j  =  1 _ C>(0’  the  matrices  defining 

the  structure  of  A.  (2  11) 

•  k index,  defined  a> 

kindex  =  ;  O]  Oj  -  o,  ].  (".(>) 

wliere  the  o{i)  are  defiieMl  b\  (2  12) 

•  indexkc  and  indexkr.  ch  liiu'd  a> 

indexkc  ^  [  c;  c-j  c,  ].  (7.7) 

indexkr  ^  !  r,  r2  ^  ]•  (”-^) 

where  c,  and  r,  are  define<l  by  (2.12). 

•  ctype,  which  defuu*s  the  performance  criterion  which  will  be 
used 


ctype 

performance  Criterion 

1 

-1 

( \»ntinuous-time  'Hn 
scaled  Popov 
Discrete-Time 

•  kl,  k2,  . .  .kr  ,  initial  values  for  the  parameter  matrices  to  be 
optimized  (2.8) 


CHAPTER  7.  PROGRAM  DESCRIPTIONS 


49 


•  N,  defined  as  N=  v 

In  addition,  scaled  Popov  Synthesis  requires  the  follow^ing  variables 
also  be  defined: 

•  Bd,Dyd,Ce,Deu,Ded,Dew,Dzd,  the  matrices  of  the  decentral¬ 
ized  static  output  feedback  realization  (2.2)  defining  the  model 
uncertainty. 

•  WW  and  ZZ,  initial  values  for  the  scaling  and  stability  multiplier 
matrices  (3.12), 

•  blk,  a  matrix  which  defines  the  block  structure  of  the  uncer¬ 
tainty  matrix  A. 

The  optimized  parameters  of  the  controller  are  written  to  the  file 
outfile. mat. 

The  optional  arguments  of  optgain  are  defineds  as  follows: 

•  outon  >  1  is  an  integer  which  determines  how  often  during  the 
iteration  process  output  will  be  displayed  to  the  screen  (and 
saved  to  the  diagnostic  files,  if  chosen;  see  below).  Screen  dis¬ 
plays  are  updated  every  outon  iterates.  The  default  value  for 
outon  is  1. 

•  diag  =  0/1  is  a  flag  for  turning  the  diagnostics  7'ecorder  on  or 
off.  The  diagnostic  recorder  saves  such  information  as  gradient 
norm  vs.  iteration  no.,  cost  function  vs.  iteration  no.,  etc. 
If  diag  =  1,  the  diagno^llc  recorder  is  on,  and  optgain  will 
generate  .mat  files  with  tins  <lata.  if  diag  =  0,  the  recorder  is 
off.  The  default  valu»'  for  diag  i>  0 

•  max  it  >  1  is  an  integer  which  represents  the  maximum  allow- 
al)le  number  of  iterations  I  he  default  value  for  maxit  is  500. 

•  method  is  a  flag  variable  which  determines  which  search  algo¬ 

rithm  the  quasi-Newton  op! imizat  ion  code  will  use  to  search  for 
the  next  candidate  parameter  \eetor.  as  follows:  (lor  more 


method 

Se  arch  .Algorithm 

1 

Line  Search 

2 

Double  Dog- Leg  Search 

3 

'‘Hook"  Step 

information  on  these  methods,  see  [4]).  The  default  value  for 
method  is  1. 

♦  tolfaO  1  is  a  multiplication  factor  for  the  code-supplied  tol¬ 
erances  within  the  cpiasi-Newton  algorithm,  which  are  on  the 
order  of  the  machine  epsilon.  By  increasing  tolfac,  the  user 
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can  relax  the  tolerances  for  the  optimization  stopping  criteria. 
The  default  value  for  tolfac  is  1. 

optgain  returns  as  output  fmin,  the  minimum  value  of  the  cost 
function  it  achieves,  the  termination  code  from  the  optimization 
routine,  info,  and  noits  the  total  number  of  iterations  performed, 
info  is  interpereted  as  follows, 


info 

Reason  for  Termination 

0 

Optimal  Solution  Found 

1 

Small  Gradient 

2 

Small  Step  Length 

3 

Unable  to  Find  Lower  Cost 

4 

Iteration  Limit  maxit  Exceeded 

5 

Too  Many  Large  Steps 
(Unbounded  Cost  Function) 

-1 

Insufficient  Workspace 

I-^xam  pies 

we  set  up  a  reduced-order,  continous-time  T^o-optimal  control  problem, 
using  a  centralized,  strictly-proper  dynamic  compensator. 

>>  A  *  [zeros(5 , 1) ,eye(5) ;-l  -2  -24  -12  -24  -4]; 

>>  B  *  [zeros  (5 , 1 ) ;  1] ;  ^ 

>>  C  e  [1  0  0  0  0  0] : 

>>  D  «  0; 

>>  D1  *  [B,zeros(6. 1)] ; 

»  D2  -  [0.1] ; 

>>  El  «  [C;zeros(  1  .€)] ; 

>>  E2  -  [0; 1] ; 

>>  EO  ■  zeros (2 , 2) ; 

»  dsoffornatC  »dataM.[l.l] , A ,B .01  .C,E1  .D.D2 ,E2,E0) ; 

»  [kl .k2.k3,k4,te8t]  •  rolqg(A.B,Dl,C,El,D,D2,E2,E0, 1) ; 

»  clear  k4 
>>  save  init  kl  k2  k3 
>>  clear 
>>  ctype  *  1; 

>>  load  data 
>>  load  init 
>>  who 
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Your  variables 

are: 

A 

Cz 

Dzv 

QL31 

indexkc 

k3 

Bu 

Dyu 

N 

QRll 

indexkr 

k index 

Bv 

Dyw 

QLll 

QR21 

kl 

Cy 

Dzu 

QL21 

QR31 

k2 

» 

save  runl 

» 

clear 

» 

[fmin,  info, 

noits]  = 

optgainC ’runl ’ 

, ’outf ile 

’) 

USING  LINE  SEARCH  ALGORITHM 


ITERATE 

FUNCTION  VALUE 

0 

0.4106206586956E+00 

1 

0.4055239785049E+00 

2 

0.35535 104 229 lOE+OO 

3 

0.3550592319726E+00 

4 

0. 35386650205 19E+00 

5 

0 . 3532465844085E+00 

24 

0.3157312171398E+00 

25 

0.3153669369866E+00 

26 

0.31S31210n808E+00 

27 

0 . 3 1 53090709009E+00 

28 

0.3153090381912E+00 

TOTAL  ITERATIONS 

00 

CM 

N 

UNCMND  WARNING  — 

•  INFO  -  1:  PROBABLY  CONVERGED,  GRADIENT  SMALL 

fmxn  « 

0.3153 
info  » 

1 

noits  ■ 

28 


» 
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pv_lmi 


Purpose 

To  provide  initializing  stability  multiplier  and  scaling  matrices  for 
scaled  Popov  synthesis. 

Synopsis 

CWW,ZZ,fail]  =  pv_lmi(A,B,C,R, gamma, blk) 

Description 

pv_lmi  finds  a  stability  multiplier  matrix  ZZ  and  scaling  matrix  WW 
for  use  with  continuous-time  scaled  Popov  criterion  performance  op- 
timization  problems.  These  matrices  correspond  to  the  solution  of 
the  scaled  Popov  Riccati  equation  (3.13)  as  the  solution  to  an  opti¬ 
mization  problem  subject  to  several  linear  matrix  inequality  (LMI) 
constraints  [2].  The  optimization  problem  is  defined  as 

min  Ir  P 

subject  to 

■  {A  +  ')-'^BC)'^P+ P{A  +  t~' RC)  + R  PB +  C'^Z +  {A  +  i-^BC)'^C'^W 
B''  P  +  ZC  +WC( A  + 'I-'  BC)  WCB  +  B'^C'^W  -~)Z 

P>D 
Z  >Q 

fail  is  a  flag  which  is  returned  0  if  pv_lmi  was  successful  in  finding 
a  feasible  pair  of  matrices,  or  1  otlierwise.  ZZ  and  WW  are  returned 
/t's  Null  if  the  LMI  solver  is  unsuc(’(*ssful. 
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rlqgly,  rlqginf,  rlqgover 


Purpose 

Find  the  reduced-order  LQG  compensator  with  an  Tioo  bound. 
Synopsis 

[Ac,  Be,  Cc,  cost]  «  rlqgly(A,B,C,D,  nc,  gammaO,  gamma,  beta,  El, 
E2,  Eli,  E2i,  Dl,  D2) 

[Ac,  Be,  Cc,  cost]  =  rlqginf  (A, B,C,D,  nc,  gammaO,  gamma,  beta,  El, 
E2,  Eli,  E2i,  Dl,  D2) 

[Ac,  Be,  Cc,  cost]  =  rlqgover(A,B,C,D,  nc,  gammaO,  gamma,  beta, 

El,  E2,  Eli,  E2i  ,  Dl,  D2) 


Desr ript  ion 

For  a  given  order  linear  })lant  with  open-loop  state  space  real¬ 
ization  given  by 


i{t) 

=  Arit)  liu(t)  -h  Diw{t) 

(7.9) 

y{i) 

zz  (>(/)  -f  Ihu  {t) 

(7.10) 

;(0 

=  K\X(i)  +  Enu(t) 

(7.11) 

c(0 

zz  xit)  A  u(t) 

(7.12) 

rlqgly.  rlqginf.  and  rlqgover  .F- .  niul  O,  a  state 

realization  for  tlie  I>Q(  J  i  Hj  optmiah  nc^^  ord<*r  roinpeusator 
vvlneh  \  ields  a  closed- loof^  with  H  ^  iionu  bounded  by  gamma. 

I  li»-  clo>e<j-loop  cost  IS  given  by  cost  gammaO  {^u)  is  the  initial  7 
arel  should  always  b<*  greater  than  gamma  beta  (.y  0)  is  a  positive 

number  'I  he  n'sulting  ccunpen>ator  from  rlqginf  is  in  the  input 
n-rmal  Hiecati  form  while  that  fr<.m  rlqgly  i>  m  by's  form. 

I  \  m  pies 


>>  A-2eroB(8); 

>>  B-2ero8(8,l) : 

>>  C«2ero8(l ,8) ; 

»  D-0; 

»  A(l:8,l)  *  [-0.161;  -6 . 004 ; -0 . 5822 ;  -9.9835;  -0.4073; 
0;  0;]; 


-3.982; 


>>  for  i  “I:?  A(i.i**‘l)  ■  1;  end 

»  B(l;8,l)  =  [0;  0;  0.0064;  0.00235;  0.0713;  1 . 0002 ; 0 . 1045 ;  0.9955;] 

>>  C(l,l)  =  1.0; 
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»  El  =  zeros(2,8) ; 

»  Eld, 1:8)  =  0.001  ♦  [0  0  0  0  0.55  11  1.32  18]; 

»  Eli  =  El; 

»  E2  =  [0;1]; 

»  E2i  =  [0;0] ; 

»  D1  -  zeros(8,2) ; 

»  Dl(l:8,l)  =  B; 

»  D2  =  [0  1]  ; 

»  [Ac, Be, Cc, cost]  =  rlqgly(A,B,C.D,2,1.0e3,3.8,100,El,E2,Eli,E2i,Dl,D2) 

Ac  = 

0  1.0000 
-0.0965  -0.2452 

Be  = 

-0.1968 
-0.1410 
Cc  = 

1.0000  -0.0000 

cost  « 

2.8821 

Algorith  m 

The  algorithms  for  rlqgly.  rlqginf.  and  rlqgover  arc  dc.scribed 
III  Chapters  I  I.  15.  Hi.  17.  and  d  ,,f  iti; 

Sec  .Also 


flqgly.  flqginf .  flqgover 
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siso2tito,  tito2siso 


Purpose 

transform  between  single-  and  dual-vector  input/output  formats 
Synopsis 

[A,Bu,Bw,Cy,Cz,Dyu,Dyw,Dzu,Dzw]  =  siso2tito(A,B,C,D,nu,ny) 
[A,B,C,D]  =  tito2siso(A,Bu,Bw,Cy,Cz,Dyu,Dyw,Dzu,Dzw) 

Description 

siso2tito  and  tito2siso  transform  continuous-  and  discrete-time 
state-space  realizations  between  single- vector  input/single- vector  out¬ 
put  (SVISVO)  and  two- vector  input /two- vector  output  (TVITVO) 
format.  Given  the  TVITVO  plant 

i:(0  =  +  + 

y(t)  =  Cyx{t)  +  I)yuU(t)  -f  Dy^.w(t), 

z{t)  =  C:x(()  +  Ihuu(t)  +  D,uw{t).  (7.13) 

where  Dyu  €  tito2siso  concatenates  the  signals  u  and  w 

to  return  the  SVISN'O  plant 

jr{t)  =  . Ij”!/ )  -r  /^u(/) 

ti{t)  -  ('xit)  ^  DuU).  (7.14) 

wUrrc  ei  ^  u  ^  ] .  and 


H  =  [ 

Hu 

IK  j 

C  = 

’• 

D  = 

■  Dyu 

Dyu  ■ 

.  Ihu 

D..  . 

Conversely,  given  the  system  (7.14).  siso2tito  will  return  the  TVITVO 
plant  (7.13)  by  selecting  the  first  n„  columns  of  B  as  the  matrix  Bu, 
the  first  iiy  rows  of  C  as  Cy,  and  l)reaking  up  D  accordingly. 

Examples 
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»  A  =  [zeros(5, 1)  ,eye(5)  ;-l  -2  -24  -12  -24  -4]; 

»  B  *  [zeros(5, 1) ;  1]  ; 

»  C  =  [1  0  0  0  0  0]  ; 

»  D  =  0; 

»  D1  =  [B,zeros(6,l)3 ; 

»  D2  =  [0,1]  ; 

»  El  =  [C;zeros(l,6)] ; 

»  E2  =  [0;  1]  ; 

»  EO  =  zeros (2, 2) ; 

»  [Abar ,Bbcu:,CbcLr,Dbar]  =  tito2siso(A,B,Dl ,C,E1  ,D,D2 ,E2 ,E0) 

Abar  = 

0  1  0  0  0  0 

0  0  1  0  0  0 

0  0  0  1  0  0 

0  0  0  0  1  0 

0  0  0  0  0  1 

-1  -2  -24  -12  -24  -4 

Bbar  * 

0  0  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 

1  1  0 

Cbar  * 

1  0  0  0  0  0 

1  0  0  0  0  0 

0  0  0  0  0  0 

Dbar  ■ 

0  0  1 

0  0  0 

1  0  0 

»  [A,B,D1,C,E1,D,D2,E2,E0]  «  s iso2ti to (Abar, Bbar , Cbar, Dbar, 1 , 1) 

A  - 

0  1  0  0  0  0 

0  0  1  0  0  0 

0  0  0  1  0  0 
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0  0  0  0  1  0 

0  0  0  0  0  1 

-1  -2  -24  -12  -24  -4 

B  = 

0 

0 

0 

0 

0 

1 

D1  = 

0  0 

0  0 

0  0 

0  0 

0  0 

1  0 

c  = 

1  0  0  0  0  0 

El  = 

1  0  0  0  0  0 

0  0  0  0  0  0 

D  « 

0 

D2  - 

0  1 

E2  - 

0 

1 

EO  - 

0  0 


0 


0 
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